#include <stdio.h>
#include <iostream>
#include <octave/oct.h>
#include <octave/file-ops.h>
#include <octave/pt-plot.h>
#include <octave/procstream.h>


static std::string Vgnuplot_command_plot;
static std::string Vgnuplot_command_replot;

static oprocstream *plot_stream = 0;

static int
set_string_var (std::string& var, const char *nm)
{
  int retval = 0;

  std::string s = builtin_string_variable (nm);

  if (s.empty ())
    {
      gripe_invalid_value_specified (nm);
      retval = -1;
    }
  else
    var = s;

  return retval;
}

static int
gnuplot_command_plot (void)
{
  return set_string_var (Vgnuplot_command_plot, "gnuplot_command_plot");
}

static int
gnuplot_command_replot (void)
{
  return set_string_var (Vgnuplot_command_replot, "gnuplot_command_replot");
}


static double mmin( const Matrix &m ) {
  double d;
  int i,j,nr,nc;

  nr = m.rows();
  nc = m.columns();

  if( nr==0 || nc==0 ) return 0.0;

  d = m(0,0);

  for(i = 0; i < nr; i++ ) {
    for( j = 0; j < nc; j++ ) {
      d = (m(i,j)<d) ? m(i,j) : d;
    }
  }
  return d;
}

static double mmax( const Matrix &m ) {
  double d;
  int i,j,nr,nc;

  nr = m.rows();
  nc = m.columns();

  if( nr==0 || nc==0 ) return 0.0;

  d = m(0,0);

  for(i = 0; i < nr; i++ ) {
    for( j = 0; j < nc; j++ ) {
      d = (m(i,j)>d) ? m(i,j) : d;
    }
  }
  return d;
}


DEFUN_DLD( __cntr__, args, , 
    "__cntr__ (x,y,z,levels): plot contour lines of a function\n\n"
    "Plot contour lines by interpolating the values of a function\n"
    "on a grid.  If x,y,z are matrices with the same dimensions\n"
    "the the corresponding elements represent points on the function.\n"
    "levels contains a vector of heights at which to draw level curves.\n" )
{
  octave_value_list retval;
  int nargin = args.length();
  int nr, nc, nl, i, j, k;
  int flag=0;
  double x11, x21, x12, x22,
  y11, y21, y12, y22,
  z11, z21, z12, z22;
  double c, x, y;
  double minx, miny, maxx, maxy;
  std::string fname;
  FILE *f;
  char cmd[256];

  if( nargin != 4 ) {
    print_usage( "cnt" );
    return retval;
  }

  // Initialize these strings
  gnuplot_command_plot();
  gnuplot_command_replot();

  Matrix xx = args(0).matrix_value();
  Matrix yy = args(1).matrix_value();
  Matrix zz = args(2).matrix_value();
  ColumnVector level = ColumnVector( args(3).vector_value());

  nl = level.length();
  nr = xx.rows();
  nc = xx.columns();

  minx = mmin(xx);
  maxx = mmax(xx);
  miny = mmin(yy);
  maxy = mmax(yy);

  if( plot_stream && *plot_stream ) {
    delete plot_stream;
    plot_stream = NULL;
  }

  plot_stream = new oprocstream( "gnuplot" );

  if( !plot_stream || ! *plot_stream )  {
    error( "Could not open pipe." );
    return retval;
  }
  
  for( k = 0; k < nl; k++ ) {
    c = level(k);
    fname = file_ops::tempnam ("", "oct-");
    f = fopen( fname.c_str(), "w" );
    mark_for_deletion( fname );

    for( i = 0; i < nr-1; i++ ) {
      for( j = 0; j < nc-1; j++ ) {
	x11 = xx(i,j );
	x21 = xx(i,j+1);
	x12 = xx(i+1,j);
	x22 = xx(i+1,j+1);

	y11 = yy(i,j );
	y21 = yy(i,j+1);
	y12 = yy(i+1,j);
	y22 = yy(i+1,j+1);

	z11 = zz(i,j );
	z21 = zz(i,j+1);
	z12 = zz(i+1,j);
	z22 = zz(i+1,j+1);

	// find bottom intersection with contour
	if(( z11>=c && z21<c ) || ( z11<c && z21>=c )) {
	  x = (c - z11)*(x21-x11)/(z21-z11) + x11;
	  y = (c - z11)*(y21-y11)/(z21-z11) + y11;
	  fprintf( f, "%lf %lf\n", x, y );
	  flag = 1;
	}

	// find top intersection with contour line
	if(( z12>=c && z22<c ) || ( z12<c && z22>=c )) {
	  x = (c - z12)*(x22-x12)/(z22-z12) + x12;
	  y = (c - z12)*(y22-y12)/(z22-z12) + y12;
	  fprintf( f, "%f %f\n", x, y );
	  flag = 1;
	}

	// find left intersection with contour line
	if(( z11>=c && z12<c ) || ( z11<c && z12>=c )) {
	  x = (c - z11)*(x12-x11)/(z12-z11) + x11;
	  y = (c - z11)*(y12-y11)/(z12-z11) + y11;
	  fprintf( f, "%f %f\n", x, y );
	  flag = 1;
	}

	// find right intersection with contour line
	if(( z21>=c && z22<c ) || ( z21<c && z22>=c )) {
	  x = (c - z21)*(x22-x21)/(z22-z21) + x21;
	  y = (c - z21)*(y22-y21)/(z22-z21) + y21;
	  fprintf( f, "%f %f\n", x, y );
	  flag = 1;
	}

	if( flag == 1 ) {
	  fprintf( f, "\n" );
	  flag = 0;
	}
      }
    }
    fclose( f );
    if( k == 0 ) {
      snprintf( cmd, sizeof(cmd), "%s [%lf:%lf] [%lf:%lf] \"%s\" title \"%g\" with lines 1\n",
	  Vgnuplot_command_plot.c_str(),
	  minx, maxx, miny, maxy,
	  fname.c_str(), c );
      *plot_stream << cmd;
    } else {
      snprintf( cmd, sizeof(cmd), "%s \"%s\" title \"%g\" with lines %d\n",
	  Vgnuplot_command_replot.c_str(),
	  fname.c_str(), c, k+1 );
      *plot_stream << cmd;
    }
  }

  plot_stream->flush();

  return retval;
}