## Copyright (C) 1999 Peter Hopfgartner
##
## This file is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## This file is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, write to the Free
## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
## 02111-1307, USA.

## Author: Peter Hopfgartner <hopfgartner@rolmail.net>
## Created: November 1999

function superficie (XX, YY, ZZ)

  if (nargin == 1)
    ## Plot on uniform grid
    ## assume that the only argument is the elevation
    ZZ = XX;
    xrow = rows (ZZ);
    xcol = columns (ZZ);
    if (is_matrix(ZZ))
      ## Set Up superficie
      ## Give plotting instruction
      tmp_in_file = tmpnam()
      [fid, msg] = fopen(tmp_in_file, "w");
      if (fid == -1)
	error("while opening superficie input file %s: \n\t%s", \
	      tmp_in_file, msg);
      endif

      fprintf ( fid, "\nQUADS");
      fprintf ( fid, "\nvertex");      
      fprintf ( fid, "\n%d", 3 * xrow * xcol);
      fprintf ( fid, "\n  {");
      for i = 1:xrow
	for j = 1:xcol
	  fprintf ( fid, "\n     %f, %f, %f", i - 1, j - 1, XX(i, j));
	  if (i != xrow || j != xcol)
	    fprintf ( fid, ",");
	  endif
	endfor
      endfor
      fprintf ( fid, "\n  }");

      fprintf ( fid, "\nindex");
      fprintf ( fid, "\n%d", 4 * (xrow-1) * (xcol-1));
      fprintf ( fid, "\n  {");
      for i=1:xrow-1
	for j=1:xcol-1
	  i1 = (i - 1) * xcol + j - 1;
	  i2 = (i - 1) * xcol + j;
	  i3 = i * xcol + j;
	  i4 = i * xcol + j - 1;
	  fprintf ( fid, "\n      %d, %d, %d, %d", i1, i2, i3, i4);
	  if (i != xrow - 1 || j != xcol - 1)
	    fprintf ( fid, ",");
	  endif
	endfor
      endfor
      fprintf (fid, "\n  }");
      fflush (fid);
      fclose(fid);
      command_string = sprintf('superficie  %s > superficie.log',tmp_in_file)
      system (command_string, 'sync');
      command_string = sprintf("rm %s", tmp_in_file);
      system(command_string);
    else
      error ("argument must be a matrix");
    endif

  elseif (nargin == 2)
    ## Plot mesh only, assume that the two arguments are the coordinates
    ## of the mesh
    xrow = rows (XX);
    xcol = columns (XX);
    if (xrow == rows (YY) && xcol == columns (YY))
      tmp_in_file = tmpnam()
      [fid, msg] = fopen(tmp_in_file, "w");
      if (fid == -1)
	error("while opening spline input file in spline_int: \n\t%s", msg);
      endif

      fprintf ( fid, "\nQUADS");
      fprintf ( fid, "\nvertex");      
      fprintf ( fid, "\n%d", 3 * xrow * xcol);
      fprintf ( fid, "\n  {");
      for i = 1:xrow
	for j = 1:xcol
	  fprintf ( fid, "\n     %f, %f, %f", XX(i, j), YY(i, j), 0.);
	  if (i != xrow || j != xcol)
	    fprintf ( fid, ",");
	  endif
	endfor
      endfor
      fprintf ( fid, "\n  }");

      fprintf ( fid, "\nindex");
      fprintf ( fid, "\n%d", 4 * (xrow-1) * (xcol-1));
      fprintf ( fid, "\n  {");
      for i=1:xrow-1
	for j=1:xcol-1
	  i1 = (i - 1) * xcol + j - 1;
	  i2 = (i - 1) * xcol + j;
	  i3 = i * xcol + j;
	  i4 = i * xcol + j - 1;
	  fprintf ( fid, "\n      %d, %d, %d, %d", i1, i2, i3, i4);
	  if (i != xrow - 1 || j != xcol - 1)
	    fprintf ( fid, ",");
	  endif
	endfor
      endfor
      fprintf (fid, "\n  }");
      fflush (fid);
      fclose(fid);
      command_string = sprintf('superficie  %s > superficie.log',tmp_in_file)
      system (command_string, 'sync');
      command_string = sprintf("rm %s", tmp_in_file);
      system(command_string);
    else
      error ("dimension of matrices must match");
    endif
  elseif (nargin == 3)
    ## Plot mesh and elevation, assume that the first two arguments
    ## are the coordinates of the mesh and the third is the elevation
    xrow = rows (XX);
    xcol = columns (XX);
    if (xrow == rows (YY) && xcol == columns (YY) &&\
	xrow == rows (ZZ) && xcol == columns (ZZ))

      tmp_in_file = tmpnam()
      [fid, msg] = fopen(tmp_in_file, "w");
      if (fid == -1)
	error("while opening spline input file in spline_int: \n\t%s", msg);
      endif

      fprintf ( fid, "\nQUADS");
      fprintf ( fid, "\nvertex");      
      fprintf ( fid, "\n%d", 3 * xrow * xcol);
      fprintf ( fid, "\n  {");
      for i = 1:xrow
	for j = 1:xcol
	  fprintf ( fid, "\n     %f, %f, %f", XX(i, j), YY(i, j), ZZ(i, j));
	  if (i != xrow || j != xcol)
	    fprintf ( fid, ",");
	  endif
	endfor
      endfor
      fprintf ( fid, "\n  }");

      fprintf ( fid, "\nindex");
      fprintf ( fid, "\n%d", 4 * (xrow-1) * (xcol-1));
      fprintf ( fid, "\n  {");
      for i=1:xrow-1
	for j=1:xcol-1
	  i1 = (i - 1) * xcol + j - 1;
	  i2 = (i - 1) * xcol + j;
	  i3 = i * xcol + j;
	  i4 = i * xcol + j - 1;
	  fprintf ( fid, "\n      %d, %d, %d, %d", i1, i2, i3, i4);
	  if (i != xrow - 1 || j != xcol - 1)
	    fprintf ( fid, ",");
	  endif
	endfor
      endfor
      fprintf (fid, "\n  }");
      fflush (fid);
      fclose(fid);
      command_string = sprintf('superficie  %s > superficie.log',tmp_in_file)
      system (command_string, 'sync');
      command_string = sprintf("rm %s", tmp_in_file);
      system(command_string);
    else
      error ("dimension of matrices must match");
    endif
  else
    error ("1, 2 or 3 arguments required");
  endif
endfunction