--- /dev/null
+## Copyright (C) 2003-2012 Shai Ayal
+##
+## This file is part of Octave.
+##
+## Octave 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 3 of the License, or (at
+## your option) any later version.
+##
+## Octave 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, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{c}, @var{lev}] =} contourc (@var{x}, @var{y}, @var{z}, @var{vn})
+## Compute isolines (contour lines) of the matrix @var{z}.
+## Parameters @var{x}, @var{y} and @var{vn} are optional.
+##
+## The return value @var{lev} is a vector of the contour levels.
+## The return value @var{c} is a 2 by @var{n} matrix containing the
+## contour lines in the following format
+##
+## @example
+## @group
+## @var{c} = [lev1, x1, x2, @dots{}, levn, x1, x2, ...
+## len1, y1, y2, @dots{}, lenn, y1, y2, @dots{}]
+## @end group
+## @end example
+##
+## @noindent
+## in which contour line @var{n} has a level (height) of @var{levn} and
+## length of @var{lenn}.
+##
+## If @var{x} and @var{y} are omitted they are taken as the row/column
+## index of @var{z}. @var{vn} is either a scalar denoting the number of lines
+## to compute or a vector containing the values of the lines. If only one
+## value is wanted, set @code{@var{vn} = [val, val]};
+## If @var{vn} is omitted it defaults to 10.
+##
+## For example:
+##
+## @example
+## @group
+## x = 0:2;
+## y = x;
+## z = x' * y;
+## contourc (x, y, z, 2:3)
+## @result{} 2.0000 2.0000 1.0000 3.0000 1.5000 2.0000
+## 2.0000 1.0000 2.0000 2.0000 2.0000 1.5000
+## @end group
+## @end example
+## @seealso{contour}
+## @end deftypefn
+
+## Author: Shai Ayal <shaiay@users.sourceforge.net>
+
+function [cout, lev] = contourc (varargin)
+
+ if (nargin == 1)
+ vn = 10;
+ z = varargin{1};
+ [nr, nc] = size (z);
+ x = 1:nc;
+ y = 1:nr;
+ elseif (nargin == 2)
+ vn = varargin{2};
+ z = varargin{1};
+ [nr, nc] = size (z);
+ x = 1:nc;
+ y = 1:nr;
+ elseif (nargin == 3)
+ vn = 10;
+ x = varargin{1};
+ y = varargin{2};
+ z = varargin{3};
+ elseif (nargin == 4)
+ vn = varargin{4};
+ x = varargin{1};
+ y = varargin{2};
+ z = varargin{3};
+ else
+ print_usage ();
+ endif
+
+ if (!ismatrix (z) || isvector (z) || isscalar (z))
+ error ("contourc: Z argument must be a matrix");
+ endif
+
+ if (isscalar (vn))
+ vv = linspace (min (z(:)), max (z(:)), vn+2)(2:end-1);
+ else
+ vv = unique (sort (vn));
+ endif
+
+ if (isvector (x) && isvector (y))
+ c = __contourc__ (x(:)', y(:)', z, vv);
+ else
+ ## Indexes x,y for the purpose of __contourc__.
+ ii = 1:size (z,2);
+ jj = 1:size (z,1);
+
+ ## Now call __contourc__ for the real work...
+ c = __contourc__ (ii, jj, z, vv);
+
+ ## Map the contour lines from index space (i,j) back
+ ## to the original grid (x,y)
+ i = 1;
+
+ while (i < size (c,2))
+ clen = c(2, i);
+ ind = i + [1 : clen];
+
+ ci = c(1, ind);
+ cj = c(2,ind);
+
+ ## due to rounding errors some elements of ci and cj
+ ## can fall out of the range of ii and jj and interp2 would
+ ## return NA for those values.
+ ## The permitted range is enforced here:
+
+ ci = max (ci, 1); ci = min (ci, size (z, 2));
+ cj = max (cj, 1); cj = min (cj, size (z, 1));
+
+ c(1, ind) = interp2 (ii, jj, x, ci, cj);
+ c(2, ind) = interp2 (ii, jj, y, ci, cj);
+
+ i = i + clen + 1;
+ endwhile
+ endif
+
+ if (nargout > 0)
+ cout = c;
+ lev = vv;
+ endif
+
+endfunction
+
+%!test
+%! x = 0:2;
+%! y = x;
+%! z = x' * y;
+%! [c_actual, lev_actual]= contourc (x, y, z, 2:3);
+%! c_expected = [2, 1, 1, 2, 2, 3, 1.5, 2; 4, 2, 2, 1, 1, 2, 2, 1.5];
+%! lev_expected = [2 3];
+%! assert (c_actual, c_expected, eps)
+%! assert (lev_actual, lev_expected, eps)
+
+