#include "sys/types.h"
#include "sys/stat.h"
#include "stdlib.h"
#include "stdio.h"
#include "math.h"
#include "unistd.h"
#include "mart_plot.h"

#define TEMPFILE "/tmp/tempfileXXXXXX"
#define PANIC(a) do {\
	           perror(a);\
		   if (temp_name) unlink(temp_name);\
                   exit(1);\
		 } while(0)


void mart_plot (double *u, double **mesh)
{
  int i, j;
  deform_grad gradU;
  deform u_local;
  double **grad, **C, **U0, **U1, **A0, **A1;
  double scale;			/* grey-scale for plotting */

  grad = dmatrix (1, 2, 1, 2);
  C = dmatrix (1, 2, 1, 2);
  U0 = dmatrix (1, 2, 1, 2);
  U1 = dmatrix (1, 2, 1, 2);
  A0 = dmatrix (1, 2, 1, 2);
  A1 = dmatrix (1, 2, 1, 2);

  for (i = 0; i < N; i++)
    for (j = 0; j < N; j++)
      {
	u_local.u0[0][0] = *(u + 2 * i * (N + 1) + 2 * j + 1);
	u_local.u1[0][0] = *(u + 2 * i * (N + 1) + 2 * j + 2);
	u_local.u0[1][0] = *(u + 2 * (i + 1) * (N + 1) + 2 * j + 1);
	u_local.u1[1][0] = *(u + 2 * (i + 1) * (N + 1) + 2 * j + 2);
	u_local.u0[0][1] = *(u + 2 * i * (N + 1) + 2 * (j + 1) + 1);
	u_local.u1[0][1] = *(u + 2 * i * (N + 1) + 2 * (j + 1) + 2);
	u_local.u0[1][1] = *(u + 2 * (i + 1) * (N + 1) + 2 * (j + 1) + 1);
	u_local.u1[1][1] = *(u + 2 * (i + 1) * (N + 1) + 2 * (j + 1) + 2);

	/* compute deformation gradients at the center of each element */
	find_gradient (&u_local, &gradU);
	grad[1][1] = gradU.grad0[0];
	grad[2][1] = gradU.grad1[0];
	grad[1][2] = gradU.grad0[1];
	grad[2][2] = gradU.grad1[1];

//      cout << "\n" << grad[1][2] << " " << grad[2][1] << " " << grad[1][1] << " " << grad[2][2];

	/* compute matrix C = FF */
	C[1][1] = DSQR (grad[1][1]) + DSQR (grad[2][1]);
	C[1][2] = grad[1][1] * grad[1][2] + grad[2][1] * grad[2][2];
	C[2][1] = grad[1][1] * grad[1][2] + grad[2][1] * grad[2][2];
	C[2][2] = DSQR (grad[1][2]) + DSQR (grad[2][2]);

//      cout << c[1][2] << " " << c[2][1] << " " << c[1][1] << " " << c[2][2] << "\n";

	/* contruct two-well */
	U0[1][1] = 1.0;
	U0[1][2] = epsilon;
	U0[2][1] = epsilon;
	U0[2][2] = 1.0 + DSQR (epsilon);
	U1[1][1] = 1.0;
	U1[1][2] = -1 * epsilon;
	U1[2][1] = -1 * epsilon;;
	U1[2][2] = 1.0 + DSQR (epsilon);

	/* compute the Grey-Scale */
	matrix_sum (0, 2, 2, C, U0, A0);	// A0 = C - U0
	matrix_sum (0, 2, 2, C, U1, A1);	// A1 = C - U1

//      cout << "\n" << matrix_norm(2,2,A0) << " " << matrix_norm(2,2,A1);
	scale =
	  matrix_norm (2, 2,
		       A0) / (matrix_norm (2, 2, A0) + matrix_norm (2, 2,
								    A1));

//      if (scale < TOL_DIST)
//        *(mesh[N - j] + (i + 1)) = 0.0;
//      else if (fabs (scale - 1) <= TOL_DIST)
//        *(mesh[N - j] + (i + 1)) = 1.0;
//      else
//        *(mesh[N - j] + (i + 1)) = .5;

	if (scale <= 1 && scale > 0.75)
	  *(mesh[N - j] + (i + 1)) = 1.0;
	else if (scale <= 0.75 && scale > 0.5)
	  *(mesh[N - j] + (i + 1)) = 0.75;
	else if (scale <= 0.5 && scale > 0.25)
	  *(mesh[N - j] + (i + 1)) = 0.5;
	else if (scale <= 0.25 && scale >= 0)
	  *(mesh[N - j] + (i + 1)) = .25;
	else
	  *(mesh[N - j] + (i + 1)) = 0.0;

      }


  /*********************************************************************/
  /* Call gnuplot in c program */

  FILE *command, *data;
  char *temp_name;
//  char temp_name[128] = TEMPFILE;
  double a, b;
  int temp_fd;

  if ((temp_name = tmpnam ((char *) 0)) == 0)
    PANIC ("tmpnam failed");
  if (mkfifo (temp_name, S_IRUSR | S_IWUSR) != 0)
    PANIC ("mkfifo failed");

//  if ((temp_fd = mkstemp(temp_name)) == -1)
//    PANIC ("mkstemp failed");

  command = popen ("gnuplot", "w");
  fprintf (command, "splot \"%s\" title \"2-D Simulation\" with lines\n",
	   temp_name);
  fflush (command);
  data = fopen (temp_name, "w");
  for (i = 1; i <= N; i++)
    {
      for (j = 1; j <= N; j++)
	{
	  fprintf (data, "%d %d %f\n", i, j, mesh[i][j]);
	}
    }
  fclose (data);
  fprintf (stderr, "press enter to continue...");
  fflush (stderr);
  getchar ();

  fclose (command);
  unlink (temp_name);

  return;
}