function contour( x, y, z, levels )
# contour (x,y,z,levels): plot contour lines of a function
#
#  Plot contour lines by interpolating the values of a function
#  on a grid.  If x,y,z are matrices with the same dimensions
#  the the corresponding elements represent points on the function.
#  If levels is a scalar, it specifies a number of contour lines
#  to draw at evenly spaced heights.  It can also be a vector of 
#  specific heights at which to draw level curves.
#
#  If x is a vector and y is a vector, then z must be a matrix
#  with the same number of rows as x and the same number of
#  columns as y.  The corresponding points from which the contour
#  lines are interpolated are x(i),y(j),z(i,y).

if( nargin < 3 )
  usage( "contour(x,y,z,levels)" );
  return;
endif
 
# If levels is unspecified, choose 5 level curves.
if( nargin == 3 )
  levels = 5;
endif

# If levels is a scalar, choose that many equally spaced curves.
if( is_scalar( levels ) )
  minz = min(min(z));
  maxz = max(max(z));
  diff = (maxz-minz)/(2*levels);

  levels = linspace(minz+diff,maxz-diff,levels );
endif

if( !is_vector( levels ))
  error( "levels must be a vector or a scalar." );
  return
endif

# If vector arguments were used, convert them to a
# rectangular mesh first.
if( is_vector(x) && is_vector(y))
  if( length(x) != rows(z) )
    error( "z must have the same number of rows as the length of x" );
    return
  endif

  if( length(y) != columns(z) )
    error( "z must have the same number of columns as the length of y" );
    return
  endif

  [x,y] = meshgrid(x,y);
endif

if( rows(x) != rows(z) || rows(y) != rows(z) ||
    columns(x) != columns(z) || columns(y) != columns(z) )
    error( "Matrix arguments must all have the same dimensions" );
    return;
endif

# Save the ishold state since we need to turn it on to
# draw multiple contour lines.
oldhold = ishold;

__cntr__(x,y,z,levels);

# reset the ishold state if it was off
if( !oldhold )
  hold off;
endif

endfunction