#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "visit_writer.h"

#define STRINGLENGTH 255
#define NDIM_READ 2
#define SET 1
#define TWO_DIM 2
#define THREE_DIM 3


int read_header(char *filename, char **varnames, int *vel_vec_dim, double **minmax);

int
main (int argc, char *argv[])
{

  FILE *fp, *fout;
  int StreamLength, StreamLengthCheck;
  char jobLabel[STRINGLENGTH];
  int numelm, npoints, nfreq, nus_sn;
  int NUM_DIM, dimension, segments;
  int ntime;
  int NumbScalarVariables, NumbVectorVariables, NumbTensorVariables;
  int nzmax;
  double time;
  int i, iii, j, ibreak, nc, nv, *nodes, **nump, inode;
  float **pData, *pptr;
  float *vptr;
  float *pDataR;
  double deg;
  double *ScalarVariables, *VarScalarMinMax;
  double *VectorVariables, *VarVectorMinMax;
  double **Variables;
  double *ScalarVariablesDummy, *VarScalarMinMaxDummy;
  double **VectorVariablesDummy, **VarVectorMinMaxDummy;
  double **minmax;
  double minmaxvalue;
  int *NodeCount;
  char **varnames;
  int nvars;
  char InString[STRINGLENGTH];
  int *conn, indx;
  int *vardims, *centering;
  int *celltypes;
  int which_scalar = 1;
  int which_vector = 0;
  float pi_factor = 2;
  int vel_vec_dim = 3;
  char outputfilename[STRINGLENGTH];


  dimension = TWO_DIM; 
  segments = 1;
  NUM_DIM = 2;

  fp = fopen (argv[1], "r");

  fread (&StreamLength, sizeof (int), 1, fp);
  fread (&InString, sizeof (char), StreamLength, fp);
  fread (&StreamLengthCheck, sizeof (int), 1, fp);

  if (StreamLength != StreamLengthCheck)
    {
      printf ("HEADER START string read incorrectly\n");
      printf ("StreamLength = %d \n", StreamLength);
      printf ("StreamLength = %d \n", StreamLengthCheck);
      exit(1);	
    }

  fread (&StreamLength, sizeof (int), 1, fp);	/* Check version */
  fread (&InString, sizeof (char), StreamLength, fp);
  fread (&StreamLengthCheck, sizeof (int), 1, fp);

  if (StreamLength != StreamLengthCheck)
    {
      printf ("VERSION read incorrectly\n");
      printf ("StreamLength = %d \n", StreamLength);
      printf ("StreamLength = %d \n", StreamLengthCheck);
      exit(1);	
    }

  /* 
   * Read the Job Label 
   */
  fread (&StreamLength, sizeof (int), 1, fp);
  fread (jobLabel, sizeof (char), StreamLength, fp);
  fread (&StreamLengthCheck, sizeof (int), 1, fp);

  if (StreamLength != StreamLengthCheck)
    {
      printf ("Filename read incorrectly\n");
      printf ("StreamLength = %d \n", StreamLength);
      printf ("StreamLength = %d \n", StreamLengthCheck);
      exit(1);	
    }

  /* 
   * Read simulation time steps: ntime
   * Number of quads: numelm
   * Number of points: npoints
   * Don't know what: nus_sn
   * Number of Frequencies: nfreq 
   */
  fread (&StreamLength, sizeof (int), 1, fp);
  fread (&ntime, sizeof (ntime), 1, fp);
  fread (&numelm, sizeof (int), 1, fp);
  fread (&npoints, sizeof (int), 1, fp);
  fread (&nfreq, sizeof (int), 1, fp);
  fread (&nus_sn, sizeof (int), 1, fp);
  fread (&StreamLengthCheck, sizeof (int), 1, fp);

  if (StreamLength != StreamLengthCheck)
    {
      printf ("Integer parameter read incorrectly\n");
      printf ("StreamLength = %d \n", StreamLength);
      printf ("StreamLength = %d \n", StreamLengthCheck);
      exit(1);	
    }

  /* 
   * Read Time: time
   */
  fread (&StreamLength, sizeof (int), 1, fp);
  fread (&time, sizeof (double), 1, fp);
  fread (&StreamLengthCheck, sizeof (int), 1, fp);

  if (StreamLength != StreamLengthCheck)
    {
      printf ("TIME read incorrectly \n");
      printf ("StreamLength = %d \n", StreamLength);
      printf ("StreamLength = %d \n", StreamLengthCheck);
      exit(1);	
    }

  /* 
   * Read Number of Scalar Variables: NumbScalarVariables
   * Number of Flux Variables: NumbVectorVariables 
   * Number of Tensor Variables: NumbTensorVariables 
   */
  fread (&StreamLength, sizeof (int), 1, fp);
  fread (&NumbScalarVariables, sizeof (int), 1, fp);
  fread (&NumbVectorVariables, sizeof (int), 1, fp);
  fread (&NumbTensorVariables, sizeof (int), 1, fp);
  fread (&StreamLengthCheck, sizeof (int), 1, fp);

  if (StreamLength != StreamLengthCheck)
    {
      printf ("Number of variables read incorrectly\n");
      printf ("StreamLength=%d\n", StreamLength);
      printf ("StreamLengthCheck=%d\n\n", StreamLengthCheck);
      exit(1);	
    }
  nvars = NumbScalarVariables + NumbVectorVariables;

  varnames = (char **) malloc (sizeof (char *) * nvars);
  vardims = (int *) malloc (sizeof (int) * nvars);
  centering = (int *) malloc (sizeof (int) * nvars);
  for (i = 0; i < NumbScalarVariables; i++)
    {
      varnames[i] = (char *) malloc (sizeof (char) * STRINGLENGTH);
      centering[i] = 1;
      vardims[i] = 1;
    }
  for (i = NumbScalarVariables; i < NumbScalarVariables+NumbVectorVariables; i++)
    {
      varnames[i] = (char *) malloc (sizeof (char) * STRINGLENGTH);
      centering[i] = 1;
      vardims[i] = 3;
    }
  minmax = (double **) malloc(sizeof(double)*nvars);
  for (i=0; i < nvars; i++){
    minmax[i] = (double *) malloc(2*sizeof(double));
  }
  read_header(argv[2], varnames, &vel_vec_dim, minmax);
  /* 
   * Read nzmax
   */
  fread (&StreamLength, sizeof (int), 1, fp);
  fread (&nzmax, sizeof (int), 1, fp);
  fread (&StreamLengthCheck, sizeof (int), 1, fp);

  if (StreamLength != StreamLengthCheck)
    {
      printf ("NZMAX incorrectly read\n");
      printf ("StreamLength=%d\n", StreamLength);
      printf ("StreamLengthCheck=%d\n\n", StreamLengthCheck);
      exit(1);	
    }


  /* 
   * Read number of corners for each cell: nodes
   */
  fread (&StreamLength, sizeof (int), 1, fp);

  nodes = (int *) calloc (numelm, sizeof (int));
  for (i = 0; i < numelm; i++)
    {
      fread (nodes + i, sizeof (int), 1, fp);
    }
  fread (&StreamLengthCheck, sizeof (int), 1, fp);

  if (StreamLength != StreamLengthCheck)
    {
      printf ("\nNODES read incorrectly\n");
      printf ("StreamLength=%d\n", StreamLength);
      printf ("StreamLengthCheck=%d\n", StreamLengthCheck);
      exit(1);	
    }

  /* 
   * Read Connections
   */

  fread (&StreamLength, sizeof (int), 1, fp);

  nump = (int **) calloc (numelm, sizeof (int *));
  conn = (int *) calloc (numelm * 4, sizeof (int *));
  indx = 0;
  for (i = 0; i < numelm; i++)
    {
      nump[i] = (int *) calloc (nodes[i], sizeof (int));
      for (j = 0; j < nodes[i]; j++)
	{
	  fread (nump[i] + j, sizeof (int), 1, fp);
	  conn[indx] = *(nump[i] + j) - 1;
	  indx += 1;
	}
    }

  fread (&StreamLengthCheck, sizeof (int), 1, fp);

  if (StreamLength != StreamLengthCheck)
    {
      printf ("\nCONNECTIONS read incorrectly\n");
      printf ("StreamLength=%d\n", StreamLength);
      printf ("StreamLengthCheck=%d\n\n", StreamLengthCheck);
      exit(1);	
    }
  /* 
   * Read co-ordinates: Xp, Yp (Same for 2D and 3D)
   */
  pData = (float **) calloc (NDIM_READ, sizeof (float *));

  for (i = 0; i < NDIM_READ; i++)
    {

      fread (&StreamLength, sizeof (int), 1, fp);

      pData[i] = (float *) calloc (npoints, sizeof (float));
      for (j = 0; j < npoints; j++)
	fread (pData[i] + j, sizeof (float), 1, fp);

      fread (&StreamLengthCheck, sizeof (int), 1, fp);
    }

  if (StreamLength != StreamLengthCheck)
    {
      printf ("\nCO-ORDINATES read incorrectly\n");
      printf ("StreamLength=%d\n", StreamLength);
      printf ("StreamLengthCheck=%d\n\n", StreamLengthCheck);
      exit(1);	
    }
  fflush (stderr);
  fflush (stdout);

  pDataR = (float *) calloc (3 * segments * npoints, sizeof (float));
  pptr = pDataR;
  inode = 0;

  if (dimension == THREE_DIM)
    {
      for (j = 0; j < npoints; j++)
	{
	  *pptr++ = pData[0][j];
	  inode = inode + 1;
	  *pptr++ = 0.;
	  inode = inode + 1;

	  *pptr++ = pData[1][j];
	  inode = inode + 1;
	}
    }
  else
    {
      for (j = 0; j < npoints; j++)
	{
	  for (i = 0; i < NDIM_READ; i++)
	    {
	      *pptr++ = pData[i][j];
	      inode = inode + 1;
	    }
	  *pptr++ = 0.;
	  inode = inode + 1;
	}
    }


  /* Prepare Array for 3D View */
#if 1
  if (dimension == THREE_DIM)
    {
      deg = (pi_factor * M_PI) / segments;

      for (j = 1; j < segments; j++)
	{
	  for (i = 0; i < npoints; i++)
	    {
	      *pptr++ = cos (j * deg) * pDataR[3 * i];
	      *pptr++ = -sin (j * deg) * pDataR[3 * i];
	      *pptr++ = pDataR[3 * i + 2];
	    }
	}

    }
#endif
  /************************************************************************************
   *                                                                                  *
   * Read and assign ScalarVariables                                                  *
   *                                                                                  *
   ************************************************************************************/

  ScalarVariables = (double *) calloc (npoints, sizeof (double));
  Variables = (double **) calloc (nvars, sizeof (double *));
  for (i=0; i < NumbScalarVariables; i++){
    Variables[i] = (double *) calloc (npoints, sizeof (double));
  }
  VarScalarMinMax = (double *) calloc (2, sizeof (double));
  NodeCount = (int *) calloc (npoints, sizeof (int));
  ScalarVariablesDummy = (double *) calloc (numelm, sizeof (double));
  VarScalarMinMaxDummy = (double *) calloc (2, sizeof (double));

  for (i = 0; i < NumbScalarVariables; i++)
    {
      fread (&StreamLength, sizeof (int), 1, fp);	/* Read Variable */
      for (j = 0; j < 2; j++)
	{
	  fread (VarScalarMinMaxDummy + j, sizeof (double), 1, fp);
	}
      for (j = 0; j < numelm; j++)
	{
	  fread (ScalarVariablesDummy + j, sizeof (double), 1, fp);
	}
      fread (&StreamLengthCheck, sizeof (int), 1, fp);

      if (StreamLength != StreamLengthCheck)
	{			/* Check correct reading */
	  printf ("\nScalar variable = %d read incorrectly\n", i);
	  printf ("StreamLength=%d\n", StreamLength);
	  printf ("StreamLengthCheck=%d\n\n", StreamLengthCheck);
          exit(1);	
	}
      which_scalar = i;
      if (i == which_scalar)
	{			/* variable selector */
          minmaxvalue = minmax[i][1] - minmax[i][0];
         minmaxvalue = 1.0;
          if (minmaxvalue == 0.0)
            minmaxvalue = 1.0;
	  for (j = 0; j < 2; j++)
	    {
	      VarScalarMinMax[j] = VarScalarMinMaxDummy[j];
	    }
	  for (j = 0; j < npoints; j++)
	    {
	      ScalarVariables[j] = 0.;
	      NodeCount[j] = 0;
	    }
	  for (j = 0; j < numelm; j++)
	    {			/* distribute cell-centered values to nodes */
	      for (nc = 0; nc < nodes[j]; nc++)
		{
		  inode = nump[j][nc] - 1;
		  ScalarVariables[inode] += ScalarVariablesDummy[j];
		  NodeCount[inode]++;
		}
	    }
	  for (j = 0; j < npoints; j++)
	    {
	      Variables[i][j] =
		(double) (ScalarVariables[j] / (NodeCount[j] * minmaxvalue));
	    }
	}
    }
  /************************************************************************************
   *                                                                                  *
   * Read and assign VectorVariables                                                  *
   *                                                                                  *
   ************************************************************************************/

  VectorVariables = (double *)calloc(3 * segments * npoints, sizeof(double));
  for (i=NumbScalarVariables; i < NumbScalarVariables+NumbVectorVariables; i++){
    Variables[i] = (double *) calloc (3 * segments *npoints, sizeof (double));
  }

  VarVectorMinMax = (double *)calloc(2, sizeof(double));
  VectorVariablesDummy = (double **)calloc(NDIM_READ, sizeof(double *));
  VarVectorMinMaxDummy = (double **)calloc(NDIM_READ, sizeof(double *));

  for(nv = 0; nv < NumbVectorVariables; nv++) {             /* Read the two components of the vector */
     minmaxvalue = minmax[nv+NumbScalarVariables][1] - minmax[NumbScalarVariables+nv][0];
     minmaxvalue = 1.0;
     if (minmaxvalue == 0.0)
        minmaxvalue = 1.0;
    ibreak = 0;
    iii = NDIM_READ;
    if(nv == 0){
      iii = vel_vec_dim;
    }
    for(i = 0; i < iii; i++) {
      if(ibreak != 0){
        break;
      }
      fread(&StreamLength, sizeof(int), 1, fp);
      VarVectorMinMaxDummy[i] = (double *)calloc(2, sizeof(double));
      for(j = 0; j < 2; j++){
        fread(VarVectorMinMaxDummy[i] + j, sizeof(double), 1, fp);
      }
      VectorVariablesDummy[i] = (double *)calloc(npoints, sizeof(double));
      for(j = 0; j < npoints; j++){
        fread(VectorVariablesDummy[i] + j, sizeof(double), 1, fp);
      }
      fread(&StreamLengthCheck, sizeof(int), 1, fp);
      if(StreamLength != StreamLengthCheck) {               /* Check correct reading */
        printf("\nVector variable = %d read incorrectly\n", i);
        printf("StreamLength=%d\n", StreamLength);
        printf("StreamLengthCheck=%d\n\n", StreamLengthCheck);
      }
    }
    which_vector = nv;
    if(nv == which_vector) {                                /* variable selector */
      for(j = 0; j < 2; j++){
        VarVectorMinMax[j] = *(VarVectorMinMaxDummy[0] + j);
      }
      inode = 0;
      indx = 0;
      if(NUM_DIM == 2) {
        for(j = 0; j < npoints; j++) {
          for(i = 0; i < NDIM_READ; i++) {
            Variables[nv+NumbScalarVariables][indx] = (double) VectorVariablesDummy[i][j]/minmaxvalue;
            indx += 1;
            inode = inode + 1;
          }
          Variables[nv+NumbScalarVariables][indx] = 0.;
          indx += 1;
          inode = inode + 1;
        }
      } 
      else {
        for(j = 0; j < npoints; j++) {
          Variables[nv+NumbScalarVariables][indx] = (double) VectorVariablesDummy[0][j]/minmaxvalue;
          indx += 1;
          inode = inode + 1;
          Variables[nv+NumbScalarVariables][indx] = 0.;
          indx += 1;
          inode = inode + 1;
          Variables[nv+NumbScalarVariables][indx] = (double) VectorVariablesDummy[1][j]/minmaxvalue;
          indx += 1;
          inode = inode + 1;
        }
#if 0
        deg = (pi_factor * M_PI) / segments;
        for(j = 1; j < segments; j++) {
          for(i = 0; i < npoints; i++) {
            Variables[nv+NumbScalarVariables][indx] = (double) (cos(j * deg) * VectorVariables[3 * i]);
            indx += 1;
            Variables[nv+NumbScalarVariables][indx] = (double) (-sin(j * deg) * VectorVariables[3 * i]);
            indx += 1;
            Variables[nv+NumbScalarVariables][indx] = (double) (VectorVariables[3 * i + 2]);
            indx += 1;
          }
        }
#endif

     } 
   }  
 } 

  celltypes = (int *) calloc (numelm, sizeof (int));
  for (i = 0; i < numelm; i++)
    {
      celltypes[i] = VISIT_QUAD;
    }
  sprintf(&outputfilename[0], "%s.vtk", argv[1]); 
  write_unstructured_mesh (outputfilename, 1, npoints, pDataR, numelm, celltypes,
			   conn, nvars, vardims, centering,
			   (const char *const *) varnames,
			   Variables);

  free (pData);
  free (pDataR);
  for (i=0; i < NumbScalarVariables+NumbVectorVariables; i++){
  	free (Variables[i]);
    free(varnames[i]);
  }
	if (Variables != NULL){
    free (Variables);
    free(varnames);
	}
  free (ScalarVariables);
  free (VarScalarMinMax);
  free (ScalarVariablesDummy);
  free (VarScalarMinMaxDummy);
  free (NodeCount);

  free(VectorVariables);
  free(VarVectorMinMax);
  free(VectorVariablesDummy);
  free(VarVectorMinMaxDummy);
  for(i = 0; i < numelm; i++){                              /* deallocate connections */
    free(nump[i]);
  }
  free(nump);
  free(nodes);


  printf ("ntime: %d\n", ntime);
  printf ("npoints: %d\n", npoints);
  printf ("time: %f\n", time);
  printf ("numelm: %d\n", numelm);
  printf
    ("NumbScalarVariables %d NumbVectorVariables %d NumbTensorVariables %d\n",
     NumbScalarVariables, NumbVectorVariables, NumbTensorVariables);
  printf ("Done reading\n");

  fclose (fp);



  return EXIT_SUCCESS;

}

#define MAXNAMELENGTH 1024 

int read_header (char *filename, char **varnames, int *vel_vec_dim, double **minmax){
	char *line = NULL;
	FILE *finfo;
	size_t len;
	ssize_t nread;
	int i, j, varnum, Number_of_Scalar_Fields_, Number_of_Vector_Fields_;
	int Number_of_Tensor_Fields_;
	char name[MAXNAMELENGTH];
        char *dataset;
	float maxvar, minvar;
  char *tokens = NULL;
  int ntok;
  char delims[] = " ";
	finfo = fopen(filename, "r");

	while ((nread = getline(&line, &len, finfo)) != -1) {

		if (strncmp(line, "Number_of_Scalar_Fields_=", strlen("Number_of_Scalar_Fields_=")) == 0){
			sscanf(line, "Number_of_Scalar_Fields_=%d", &Number_of_Scalar_Fields_);
			printf("Number_of_Scalar_Fields_= %d \n", Number_of_Scalar_Fields_);
		  for (i=0; i < Number_of_Scalar_Fields_; i++){
				if ((nread = getline(&line, &len, finfo)) != -1) {
					sscanf(line, "%d%s", &varnum, name);
          dataset = (char *) malloc(strlen(line));
          sprintf(dataset, line);
          tokens = strtok (dataset, delims);
          ntok = 0;
          while (tokens != NULL){
            tokens = strtok (NULL, delims);
            if (tokens != NULL){
              printf("tokens: %s\n", tokens);
            }
            ntok+=1;
          }
          tokens = strtok (line, delims);
          sprintf(varnames[i], "%d", i+1);
          for (j=0; j < ntok - 4; j++){
            tokens = strtok (NULL, delims);
            sprintf(varnames[i]+strlen(varnames[i]), "%s_", tokens);
          }
          tokens = strtok (NULL, delims);
          sprintf(varnames[i]+strlen(varnames[i]), "%s", tokens);
          printf("varnames[i]: %d%s\n", i, varnames[i]);

          for (j=0; j < 2; j++){
            tokens = strtok (NULL, delims);
            minmax[i][j] = atof(tokens);
          }
          printf("%20.12e %20.12e\n", minmax[i][0], minmax[i][1]);
        }
      }
		} else if (strncmp(line, "Number_of_Vector_Fields_=", strlen("Number_of_Vector_Fields_=")) == 0){
      sscanf(line, "Number_of_Vector_Fields_=%d", &Number_of_Vector_Fields_);
      printf("Number_of_Vector_Fields_= %d \n", Number_of_Vector_Fields_);
      nread = getline(&line, &len, finfo);
      sscanf(line, "Number_of_Velocity_Comp_=", &vel_vec_dim);
      for (i=Number_of_Scalar_Fields_; i < Number_of_Scalar_Fields_+Number_of_Vector_Fields_; i++){
        if ((nread = getline(&line, &len, finfo)) != -1) {
          sscanf(line, "%d%s", &varnum, name);
          dataset = (char *) malloc(strlen(line)); 
          sprintf(dataset, line);
          tokens = strtok (dataset, delims);
          ntok = 0;
          while (tokens != NULL){
            tokens = strtok (NULL, delims);
            if (tokens != NULL){
              printf("tokens: %s\n", tokens);
            } 
            ntok+=1;
          }
          tokens = strtok (line, delims);
          sprintf(varnames[i], "%d", i-Number_of_Scalar_Fields_+1);
          for (j=0; j < ntok - 4; j++){
            tokens = strtok (NULL, delims);
            sprintf(varnames[i]+strlen(varnames[i]), "%s_", tokens);
	  }
          tokens = strtok (NULL, delims);
          sprintf(varnames[i]+strlen(varnames[i]), "%s", tokens);
          printf("varnames[i]: %s\n", varnames[i]);
          
          for (j=0; j < 2; j++){
            tokens = strtok (NULL, delims);
            minmax[i][j] = atof(tokens);
	        }
          printf("%20.12e %20.12e\n", minmax[i][0], minmax[i][1]);
        } 
      }
	  } else if (strncmp(line, "Number_of_Tensor_fields_=", strlen("Number_of_Tensor_fields_=")) == 0){
      sscanf(line, "Number_of_Tensor_fields_=%d", &Number_of_Tensor_Fields_);
      printf("Number_of_Tensor_Fields_= %d \n", Number_of_Tensor_Fields_);
      for (i=Number_of_Scalar_Fields_+Number_of_Vector_Fields_; i < Number_of_Scalar_Fields_+Number_of_Vector_Fields_+Number_of_Tensor_Fields_; i++){
        if ((nread = getline(&line, &len, finfo)) != -1) {
          sscanf(line, "%d%s", &varnum, name);
          dataset = (char *) malloc(strlen(line));
          sprintf(dataset, line);
          tokens = strtok (dataset, delims);
          ntok = 0;
          while (tokens != NULL){
            tokens = strtok (NULL, delims);
            if (tokens != NULL){
              printf("tokens: %s\n", tokens);
            }
            ntok+=1;
          }
          tokens = strtok (line, delims);
          sprintf(varnames[i], "%d", i-Number_of_Scalar_Fields_+Number_of_Vector_Fields_+1);
          for (j=0; j < ntok - 4; j++){
            tokens = strtok (NULL, delims);
            sprintf(varnames[i]+strlen(varnames[i]), "%s_", tokens);
    }
          tokens = strtok (NULL, delims);
          sprintf(varnames[i]+strlen(varnames[i]), "%s", tokens);
          printf("varnames[i]: %s\n", varnames[i]);

          for (j=0; j < 2; j++){
            tokens = strtok (NULL, delims);
            minmax[i][j] = atof(tokens);
          }
          if (minmax[i][0] < 1e-34) minmax[i][0] = 0.0;
          printf("%20.12e %20.12e\n", minmax[i][0], minmax[i][1]);
			}	
     	} 
 } 
 }
  fclose(finfo);

	return EXIT_SUCCESS;

}


