--- a/NEWS	Sun Mar 30 13:02:08 2014 -0700
+++ b/NEWS	Sun Mar 30 14:18:43 2014 -0700
@@ -1,6 +1,14 @@
 Summary of important user-visible changes for version 4.2:
 ---------------------------------------------------------
 
+ ** Interpolation function changes for Matlab compatibility
+
+    The interpolation method 'cubic' is now equivalent to 'pchip'
+    for interp1, interp2, and interp3.  Previously, 'cubic' was equivalent
+    to 'spline' for interp2.  This may produce different results as 'spline'
+    has continuous 1st and 2nd derivatives while 'pchip' only has a continuous
+    1st derivative.
+
  ** Other new functions added in 4.2:
 
       dir_in_loadpath

--- a/scripts/general/bicubic.m	Sun Mar 30 13:02:08 2014 -0700
+++ b/scripts/general/bicubic.m	Sun Mar 30 14:18:43 2014 -0700
@@ -240,9 +240,9 @@
 %! z = cos (6 * xx) + sin (6 * yy);
 %! x = linspace (1, -1, 30);
 %! [xx2, yy2] = meshgrid (x);
-%! z1 = interp2 (xx, yy, z, xx2, yy2, "cubic");
+%! z1 = interp2 (xx, yy, z, xx2, yy2, "spline");
 %! z2 = interp2 (fliplr (xx), flipud (yy), fliplr (flipud(z)),
-%!               fliplr (xx2), flipud (yy2), "cubic");
+%!               fliplr (xx2), flipud (yy2), "spline");
 %! z2 = fliplr (flipud (z2));
 %! assert (z1, z2, 100 * eps ())
 

--- a/scripts/general/interp1.m	Sun Mar 30 13:02:08 2014 -0700
+++ b/scripts/general/interp1.m	Sun Mar 30 14:18:43 2014 -0700
@@ -29,28 +29,29 @@
 ## One-dimensional interpolation.
 ##
 ## Interpolate input data to determine the value of @var{yi} at the points
-## @var{xi}.  If not specified, @var{x} is taken to be the indices of @var{y}.
-## If @var{y} is a matrix or an N-dimensional array, the interpolation is
-## performed on each column of @var{y}.
+## @var{xi}.  If not specified, @var{x} is taken to be the indices of @var{y}
+## (@code{1:length (@var{y})}).  If @var{y} is a matrix or an N-dimensional
+## array, the interpolation is performed on each column of @var{y}.
 ##
-## Method is one of:
+## The interpolation @var{method} is one of:
 ##
 ## @table @asis
 ## @item @qcode{"nearest"}
-## Return the nearest neighbor
+## Return the nearest neighbor.
 ##
-## @item @qcode{"linear"}
-## Linear interpolation from nearest neighbors
+## @item @qcode{"linear"} (default)
+## Linear interpolation from nearest neighbors.
 ##
 ## @item @qcode{"pchip"}
-## Piecewise cubic Hermite interpolating polynomial
+## Piecewise cubic Hermite interpolating polynomial---shape-preserving
+## interpolation with smooth first derivative.
 ##
 ## @item @qcode{"cubic"}
-## Cubic interpolation (same as @code{pchip})
+## Cubic interpolation (same as @qcode{"pchip"}).
 ##
 ## @item @qcode{"spline"}
 ## Cubic spline interpolation---smooth first and second derivatives
-## throughout the curve
+## throughout the curve.
 ## @end table
 ##
 ## Adding '*' to the start of any method above forces @code{interp1}
@@ -61,7 +62,7 @@
 ## If @var{extrap} is the string @qcode{"extrap"}, then extrapolate values
 ## beyond the endpoints using the current @var{method}.  If @var{extrap} is a
 ## number, then replace values beyond the endpoints with that number.  When
-## unspecified, @var{extrap} defaults to NA.
+## unspecified, @var{extrap} defaults to @code{NA}.
 ##
 ## If the string argument @qcode{"pp"} is specified, then @var{xi} should not
 ## be supplied and @code{interp1} returns a piecewise polynomial object.  This 
@@ -91,16 +92,16 @@
 ## xp = [0:10];
 ## yp = sin (2*pi*xp/5);
 ## lin = interp1 (xp, yp, xf);
+## near = interp1 (xp, yp, xf, "nearest");
+## pch = interp1 (xp, yp, xf, "pchip");
 ## spl = interp1 (xp, yp, xf, "spline");
-## cub = interp1 (xp, yp, xf, "cubic");
-## near = interp1 (xp, yp, xf, "nearest");
-## plot (xf, yf, "r", xf, lin, "g", xf, spl, "b",
-##       xf, cub, "c", xf, near, "m", xp, yp, "r*");
-## legend ("original", "linear", "spline", "cubic", "nearest");
+## plot (xf,yf,"r", xf,near,"g", xf,lin,"b", xf,pch,"c", xf,spl,"m", 
+##       xp,yp,"r*");
+## legend ("original", "nearest", "linear", "pchip", "spline");
 ## @end group
 ## @end example
 ##
-## @seealso{interpft, interp2, interp3, interpn}
+## @seealso{pchip, spline, interpft, interp2, interp3, interpn}
 ## @end deftypefn
 
 ## Author: Paul Kienzle
@@ -130,17 +131,18 @@
       arg = varargin{i};
       if (ischar (arg))
         arg = tolower (arg);
-        if (strcmp ("extrap", arg))
-          extrap = "extrap";
-        elseif (strcmp ("pp", arg))
-          ispp = true;
-        elseif (strcmp (arg, "right") || strcmp (arg, "-right"))
-          rightcontinuous = true;
-        elseif (strcmp (arg, "left") || strcmp (arg, "-left"))
-          rightcontinuous = false;
-        else
-          method = arg;
-        endif
+        switch (arg)
+          case "extrap"
+            extrap = "extrap";
+          case "pp"
+            ispp = true;
+          case {"right", "-right"}
+            rightcontinuous = true;
+          case {"left", "-left"}
+            rightcontinuous = false;
+          otherwise
+            method = arg;
+        endswitch
       else
         if (firstnumeric)
           xi = arg;
@@ -177,7 +179,7 @@
 
   ## determine sizes
   if (nx < 2 || ny < 2)
-    error ("interp1: table too short");
+    error ("interp1: minimum of 2 points required in each dimension");
   endif
 
   ## check whether x is sorted; sort if not.
@@ -245,7 +247,6 @@
       endif
 
     case "linear"
-
       xx = x;
       yy = y;
       nxx = nx;
@@ -322,13 +323,13 @@
 
   endswitch
 
-  if (! ispp && ! ischar (extrap))
+  if (! ispp && isnumeric (extrap))
     ## determine which values are out of range and set them to extrap,
     ## unless extrap == "extrap".
     minx = min (x(1), x(nx));
     maxx = max (x(1), x(nx));
 
-    outliers = xi < minx | ! (xi <= maxx); # this even catches NaNs
+    outliers = xi < minx | ! (xi <= maxx);  # this even catches NaNs
     if (size_equal (outliers, yi))
       yi(outliers) = extrap;
       yi = reshape (yi, szx);
@@ -346,12 +347,13 @@
 %! clf;
 %! xf = 0:0.05:10;  yf = sin (2*pi*xf/5);
 %! xp = 0:10;       yp = sin (2*pi*xp/5);
-%! lin = interp1 (xp,yp,xf, "linear");
-%! spl = interp1 (xp,yp,xf, "spline");
-%! cub = interp1 (xp,yp,xf, "pchip");
-%! near= interp1 (xp,yp,xf, "nearest");
-%! plot (xf,yf,"r",xf,near,"g",xf,lin,"b",xf,cub,"c",xf,spl,"m",xp,yp,"r*");
-%! legend ("original", "nearest", "linear", "pchip", "spline");
+%! lin = interp1 (xp,yp,xf, 'linear');
+%! spl = interp1 (xp,yp,xf, 'spline');
+%! pch = interp1 (xp,yp,xf, 'pchip');
+%! near= interp1 (xp,yp,xf, 'nearest');
+%! plot (xf,yf,'r',xf,near,'g',xf,lin,'b',xf,pch,'c',xf,spl,'m',xp,yp,'r*');
+%! legend ('original', 'nearest', 'linear', 'pchip', 'spline');
+%! title ('Interpolation of continuous function sin (x) w/various methods');
 %! %--------------------------------------------------------
 %! % confirm that interpolated function matches the original
 
@@ -359,36 +361,50 @@
 %! clf;
 %! xf = 0:0.05:10;  yf = sin (2*pi*xf/5);
 %! xp = 0:10;       yp = sin (2*pi*xp/5);
-%! lin = interp1 (xp,yp,xf, "*linear");
-%! spl = interp1 (xp,yp,xf, "*spline");
-%! cub = interp1 (xp,yp,xf, "*cubic");
-%! near= interp1 (xp,yp,xf, "*nearest");
-%! plot (xf,yf,"r",xf,near,"g",xf,lin,"b",xf,cub,"c",xf,spl,"m",xp,yp,"r*");
-%! legend ("*original", "*nearest", "*linear", "*cubic", "*spline");
+%! lin = interp1 (xp,yp,xf, '*linear');
+%! spl = interp1 (xp,yp,xf, '*spline');
+%! pch = interp1 (xp,yp,xf, '*pchip');
+%! near= interp1 (xp,yp,xf, '*nearest');
+%! plot (xf,yf,'r',xf,near,'g',xf,lin,'b',xf,pch,'c',xf,spl,'m',xp,yp,'r*');
+%! legend ('*original', '*nearest', '*linear', '*pchip', '*spline');
+%! title ('Interpolation of continuous function sin (x) w/various *methods');
 %! %--------------------------------------------------------
 %! % confirm that interpolated function matches the original
 
 %!demo
 %! clf;
+%! fstep = @(x) x > 1;
+%! xf = 0:0.05:2;  yf = fstep (xf);
+%! xp = linspace (0,2,10);  yp = fstep (xp);
+%! pch = interp1 (xp,yp,xf, 'pchip');
+%! spl = interp1 (xp,yp,xf, 'spline');
+%! plot (xf,yf,'r',xf,pch,'b',xf,spl,'m',xp,yp,'r*');
+%! title ({'Interpolation of step function with discontinuity at x==1', ...
+%!         'Note: "pchip" is shape-preserving, "spline" (continuous 1st, 2nd derivatives) is not'});
+%! legend ('original', 'pchip', 'spline');
+
+%!demo
+%! clf;
 %! t = 0 : 0.3 : pi; dt = t(2)-t(1);
 %! n = length (t); k = 100; dti = dt*n/k;
 %! ti = t(1) + [0 : k-1]*dti;
 %! y = sin (4*t + 0.3) .* cos (3*t - 0.1);
-%! ddyc = diff (diff (interp1 (t,y,ti, "cubic")) ./dti)./dti;
-%! ddys = diff (diff (interp1 (t,y,ti, "spline"))./dti)./dti;
-%! ddyp = diff (diff (interp1 (t,y,ti, "pchip")) ./dti)./dti;
-%! plot (ti(2:end-1),ddyc,'g+', ti(2:end-1),ddys,'b*', ti(2:end-1),ddyp,'c^');
-%! legend ("cubic", "spline", "pchip");
-%! title ("Second derivative of interpolated 'sin (4*t + 0.3) .* cos (3*t - 0.1)'");
+%! ddys = diff (diff (interp1 (t,y,ti, 'spline'))./dti)./dti;
+%! ddyp = diff (diff (interp1 (t,y,ti, 'pchip')) ./dti)./dti;
+%! ddyc = diff (diff (interp1 (t,y,ti, 'cubic')) ./dti)./dti;
+%! plot (ti(2:end-1),ddys,'b*', ti(2:end-1),ddyp,'c^', ti(2:end-1),ddyc,'g+');
+%! title ({'Second derivative of interpolated "sin (4*t + 0.3) .* cos (3*t - 0.1)"', ...
+%!         'Note: "spline" has continous 2nd derivative, others do not'});
+%! legend ('spline', 'pchip', 'cubic');
 
 %!demo
 %! clf;
 %! xf = 0:0.05:10;                yf = sin (2*pi*xf/5) - (xf >= 5);
 %! xp = [0:.5:4.5,4.99,5:.5:10];  yp = sin (2*pi*xp/5) - (xp >= 5);
-%! lin = interp1 (xp,yp,xf, "linear");
-%! near= interp1 (xp,yp,xf, "nearest");
-%! plot (xf,yf,"r", xf,near,"g", xf,lin,"b", xp,yp,"r*");
-%! legend ("original", "nearest", "linear");
+%! lin = interp1 (xp,yp,xf, 'linear');
+%! near= interp1 (xp,yp,xf, 'nearest');
+%! plot (xf,yf,'r', xf,near,'g', xf,lin,'b', xp,yp,'r*');
+%! legend ('original', 'nearest', 'linear');
 %! %--------------------------------------------------------
 %! % confirm that interpolated function matches the original
 
@@ -409,7 +425,7 @@
 %! %--------------------------------------------------------
 %! % red curve is left-continuous and blue is right-continuous at x = 2
 
-##FIXME: add test for n-d arguments here
+##FIXME: add test for N-d arguments here
 
 ## For each type of interpolated test, confirm that the interpolated
 ## value at the knots match the values at the knots.  Points away
@@ -637,12 +653,12 @@
 %% Test input validation
 %!error interp1 ()
 %!error interp1 (1,2,3,4,5,6,7)
-%!error <table too short> interp1 (1,1,1, "linear")
-%!error <table too short> interp1 (1,1,1, "*nearest")
-%!error <table too short> interp1 (1,1,1, "*linear")
+%!error <minimum of 2 points required> interp1 (1,1,1, "linear")
+%!error <minimum of 2 points required> interp1 (1,1,1, "*nearest")
+%!error <minimum of 2 points required> interp1 (1,1,1, "*linear")
 %!warning <multiple discontinuities> interp1 ([1 1 1 2], [1 2 3 4], 1);
 %!error <discontinuities not supported> interp1 ([1 1],[1 2],1, "pchip")
 %!error <discontinuities not supported> interp1 ([1 1],[1 2],1, "cubic")
 %!error <discontinuities not supported> interp1 ([1 1],[1 2],1, "spline")
-%!error <invalid method> interp1 (1:2,1:2,1, "bogus")
+%!error <invalid method 'bogus'> interp1 (1:2,1:2,1, "bogus")
 

--- a/scripts/general/interp2.m	Sun Mar 30 13:02:08 2014 -0700
+++ b/scripts/general/interp2.m	Sun Mar 30 14:18:43 2014 -0700
@@ -19,63 +19,63 @@
 
 ## -*- texinfo -*-
 ## @deftypefn  {Function File} {@var{zi} =} interp2 (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi})
-## @deftypefnx {Function File} {@var{zi} =} interp2 (@var{Z}, @var{xi}, @var{yi})
-## @deftypefnx {Function File} {@var{zi} =} interp2 (@var{Z}, @var{n})
+## @deftypefnx {Function File} {@var{zi} =} interp2 (@var{z}, @var{xi}, @var{yi})
+## @deftypefnx {Function File} {@var{zi} =} interp2 (@var{z}, @var{n})
+## @deftypefnx {Function File} {@var{zi} =} interp2 (@var{z})
 ## @deftypefnx {Function File} {@var{zi} =} interp2 (@dots{}, @var{method})
-## @deftypefnx {Function File} {@var{zi} =} interp2 (@dots{}, @var{method}, @var{extrapval})
+## @deftypefnx {Function File} {@var{zi} =} interp2 (@dots{}, @var{method}, @var{extrap})
+##
+## Two-dimensional interpolation.
 ##
-## Two-dimensional interpolation.  @var{x}, @var{y} and @var{z} describe a
-## surface function.  If @var{x} and @var{y} are vectors their length
-## must correspondent to the size of @var{z}.  @var{x} and @var{y} must be
-## monotonic.  If they are matrices they must have the @code{meshgrid}
-## format.
-##
-## @table @code
-## @item interp2 (@var{x}, @var{y}, @var{Z}, @var{xi}, @var{yi}, @dots{})
-## Returns a matrix corresponding to the points described by the
-## matrices @var{xi}, @var{yi}.
+## Interpolate reference data @var{x}, @var{y}, @var{z} to determine @var{zi}
+## at the coordinates @var{xi}, @var{yi}.  The reference data @var{x}, @var{y}
+## can be matrices, as returned by @code{meshgrid}, in which case the sizes of
+## @var{x}, @var{y}, and @var{z} must be equal.  If @var{x}, @var{y} are
+## vectors describing a grid then @code{length (@var{x}) == columns (@var{z})}
+## and @code{length (@var{y}) == rows (@var{z})}.  In either case the input
+## data must be strictly monotonic.
 ##
-## If the last argument is a string, the interpolation method can
-## be specified.  The method can be @qcode{"linear"}, @qcode{"nearest"} or
-## @qcode{"cubic"}.  If it is omitted @qcode{"linear"} interpolation is
-## assumed.
+## If called without @var{x}, @var{y}, and just a single reference data matrix
+## @var{z}, the 2-D region
+## @code{@var{x} = 1:columns (@var{z}), @var{y} = 1:rows (@var{z})} is assumed.
+## This saves memory if the grid is regular and the distance between points is
+## not important.
 ##
-## @item interp2 (@var{z}, @var{xi}, @var{yi})
-## Assumes @code{@var{x} = 1:rows (@var{z})} and @code{@var{y} =
-## 1:columns (@var{z})}
+## If called with a single reference data matrix @var{z} and a refinement
+## value @var{n}, then perform interpolation over a grid where each original
+## interval has been recursively subdivided @var{n} times.  This results in
+## @code{2^@var{n}-1} additional points for every interval in the original
+## grid.  If @var{n} is omitted a value of 1 is used.  As an example, the
+## interval [0,1] with @code{@var{n}==2} results in a refined interval with
+## points at [0, 1/4, 1/2, 3/4, 1].
 ##
-## @item interp2 (@var{z}, @var{n})
-## Interleaves the matrix @var{z} n-times.  If @var{n} is omitted a value
-## of @code{@var{n} = 1} is assumed.
-## @end table
-##
-## The variable @var{method} defines the method to use for the
-## interpolation.  It can take one of the following values
+## The interpolation @var{method} is one of:
 ##
 ## @table @asis
 ## @item @qcode{"nearest"}
 ## Return the nearest neighbor.
 ##
-## @item @qcode{"linear"}
+## @item @qcode{"linear"} (default)
 ## Linear interpolation from nearest neighbors.
 ##
 ## @item @qcode{"pchip"}
-## Piecewise cubic Hermite interpolating polynomial.
+## Piecewise cubic Hermite interpolating polynomial---shape-preserving
+## interpolation with smooth first derivative.
 ##
 ## @item @qcode{"cubic"}
-## Cubic interpolation from four nearest neighbors.
+## Cubic interpolation (same as @qcode{"pchip"}).
 ##
 ## @item @qcode{"spline"}
 ## Cubic spline interpolation---smooth first and second derivatives
 ## throughout the curve.
 ## @end table
 ##
-## If a scalar value @var{extrapval} is defined as the final value, then
-## values outside the mesh as set to this value.  Note that in this case
-## @var{method} must be defined as well.  If @var{extrapval} is not
-## defined then NA is assumed.
-##
-## @seealso{interp1}
+## If @var{extrap} is the string @qcode{"extrap"}, then extrapolate values
+## beyond the endpoints using the current @var{method}.  If @var{extrap} is a
+## number, then replace values beyond the endpoints with that number.  When
+## unspecified, @var{extrap} defaults to @code{NA}.  Note that if @var{extrap}
+## is used, @var{method} must be specified as well.
+## @seealso{interp1, interp3, interpn, meshgrid}
 ## @end deftypefn
 
 ## Author:      Kai Habel <kai.habel@gmx.de>
@@ -94,10 +94,14 @@
 ##       "meshgridded") to be consistent with the help message
 ##       above and for compatibility.
 
+## FIXME: Need better input validation.
+##        E.g, interp2 (1,1,1) => A(I): index out of bounds
+
 function ZI = interp2 (varargin)
+
   Z = X = Y = XI = YI = n = [];
   method = "linear";
-  extrapval = NA;
+  extrap = NA;
 
   switch (nargin)
     case 1
@@ -120,36 +124,38 @@
       if (ischar (varargin{4}))
         [Z, XI, YI, method] = deal (varargin{:});
       else
-        [Z, n, method, extrapval] = deal (varargin{:});
+        [Z, n, method, extrap] = deal (varargin{:});
       endif
     case 5
       if (ischar (varargin{4}))
-        [Z, XI, YI, method, extrapval] = deal (varargin{:});
+        [Z, XI, YI, method, extrap] = deal (varargin{:});
       else
         [X, Y, Z, XI, YI] = deal (varargin{:});
       endif
     case 6
         [X, Y, Z, XI, YI, method] = deal (varargin{:});
     case 7
-        [X, Y, Z, XI, YI, method, extrapval] = deal (varargin{:});
+        [X, Y, Z, XI, YI, method, extrap] = deal (varargin{:});
     otherwise
       print_usage ();
   endswitch
 
   ## Type checking.
-  if (!ismatrix (Z))
+  if (! ismatrix (Z))
     error ("interp2: Z must be a matrix");
   endif
-  if (!isempty (n) && !isscalar (n))
-    error ("interp2: N must be a scalar");
+  if (! isempty (n) && ! (isscalar (n) && n >= 0 && n == fix (n)))
+    error ("interp2: N must be an integer >= 0");
   endif
-  if (!ischar (method))
+  if (! ischar (method))
     error ("interp2: METHOD must be a string");
   endif
-  if (ischar (extrapval) || strcmp (extrapval, "extrap"))
-    extrapval = [];
-  elseif (!isscalar (extrapval))
-    error ("interp2: EXTRAPVAL must be a scalar");
+  if (isnumeric (extrap) && isscalar (extrap))
+    ## Typical case
+  elseif (strcmp (extrap, "extrap"))
+    extrap = [];
+  else
+    error ('interp2: EXTRAP must be a numeric scalar or "extrap"');
   endif
 
   if (method(1) == "*")
@@ -176,9 +182,7 @@
     error ("interp2: XI, YI must be numeric");
   endif
 
-
-  if (strcmp (method, "linear") || strcmp (method, "nearest") ...
-      || strcmp (method, "pchip"))
+  if (any (strcmp (method, {"nearest", "linear", "pchip", "cubic"})))
 
     ## If X and Y vectors produce a grid from them
     if (isvector (X) && isvector (Y))
@@ -186,7 +190,7 @@
     elseif (size_equal (X, Y))
       X = X(1,:)'; Y = Y(:,1);
     else
-      error ("interp2: X and Y must be matrices of same size");
+      error ("interp2: X and Y must be matrices of equal size");
     endif
     if (columns (Z) != length (X) || rows (Z) != length (Y))
       error ("interp2: X and Y size must match the dimensions of Z");
@@ -249,10 +253,11 @@
       idx = sub2ind (size (Z), yidx+jj, xidx+ii);
       ZI = Z(idx);
 
-    elseif (strcmp (method, "pchip"))
+    elseif (strcmp (method, "pchip") || strcmp (method, "cubic"))
 
       if (length (X) < 2 || length (Y) < 2)
-        error ("interp2: pchip2 requires at least 2 points in each dimension");
+        error ("interp2: %s requires at least 2 points in each dimension",
+               method);
       endif
 
       ## first order derivatives
@@ -308,23 +313,23 @@
 
     endif
 
-    if (! isempty (extrapval))
-      ## set points outside the table to 'extrapval'
-      if (X (1) < X (end))
-        if (Y (1) < Y (end))
-          ZI (XI < X(1,1) | XI > X(end) | YI < Y(1,1) | YI > Y(end)) = ...
-                  extrapval;
+    if (! isempty (extrap))
+      ## set points outside the table to 'extrap'
+      if (X(1) < X(end))
+        if (Y(1) < Y(end))
+          ZI(XI < X(1,1) | XI > X(end) | YI < Y(1,1) | YI > Y(end)) = ...
+                  extrap;
         else
-          ZI (XI < X(1) | XI > X(end) | YI < Y(end) | YI > Y(1)) = ...
-                  extrapval;
+          ZI(XI < X(1) | XI > X(end) | YI < Y(end) | YI > Y(1)) = ...
+                  extrap;
         endif
       else
-        if (Y (1) < Y (end))
-          ZI (XI < X(end) | XI > X(1) | YI < Y(1) | YI > Y(end)) = ...
-                  extrapval;
+        if (Y(1) < Y(end))
+          ZI(XI < X(end) | XI > X(1) | YI < Y(1) | YI > Y(end)) = ...
+                  extrap;
         else
-          ZI (XI < X(1,end) | XI > X(1) | YI < Y(end) | YI > Y(1)) = ...
-                  extrapval;
+          ZI(XI < X(1,end) | XI > X(1) | YI < Y(end) | YI > Y(1)) = ...
+                  extrap;
         endif
       endif
     endif
@@ -335,14 +340,13 @@
     if (isvector (X) && isvector (Y))
       X = X(:).';
       Y = Y(:);
-      if (!isequal ([length(Y), length(X)], size(Z)))
+      if (! isequal ([length(Y), length(X)], size (Z)))
         error ("interp2: X and Y size must match the dimensions of Z");
       endif
-    elseif (!size_equal (X, Y))
+    elseif (! size_equal (X, Y))
       error ("interp2: X and Y must be matrices of equal size");
-      if (! size_equal (X, Z))
-        error ("interp2: X and Y size must match the dimensions of Z");
-      endif
+    elseif (! size_equal (X, Z))
+      error ("interp2: X and Y size must match the dimensions of Z");
     endif
 
     ## Check dimensions of XI and YI
@@ -354,75 +358,84 @@
       error ("interp2: XI and YI must be matrices of equal size");
     endif
 
+    ## FIXME: Previously used algorithm for cubic.
+    ##        This produced results within a few eps of "spline".
+    ##        Matlab compatibility requires "cubic" to be a C1 algorithm
+    ##        equivalent to "pchip" so this was commented out 2014/03/30.
+    ##        This can be removed completely in the future if no problems are
+    ##        encountered.
+    #{
     if (strcmp (method, "cubic"))
       if (isgriddata (XI) && isgriddata (YI'))
-        ZI = bicubic (X, Y, Z, XI (1, :), YI (:, 1), extrapval);
+        ZI = bicubic (X, Y, Z, XI (1, :), YI (:, 1), extrap);
       elseif (isgriddata (X) && isgriddata (Y'))
         ## Allocate output
         ZI = zeros (size (X));
 
         ## Find inliers
-        inside = !(XI < X (1) | XI > X (end) | YI < Y (1) | YI > Y (end));
+        inside = !(XI < X(1) | XI > X(end) | YI < Y(1) | YI > Y(end));
 
         ## Scale XI and YI to match indices of Z
-        XI = (columns (Z) - 1) * (XI - X (1)) / (X (end) - X (1)) + 1;
-        YI = (rows (Z) - 1) * (YI - Y (1)) / (Y (end) - Y (1)) + 1;
+        XI = (columns (Z) - 1) * (XI - X(1)) / (X(end) - X(1)) + 1;
+        YI = (rows (Z) - 1) * (YI - Y(1)) / (Y(end) - Y(1)) + 1;
 
         ## Start the real work
         K = floor (XI);
         L = floor (YI);
 
         ## Coefficients
-        AY1  = bc ((YI - L + 1));
-        AX1  = bc ((XI - K + 1));
-        AY0  = bc ((YI - L + 0));
-        AX0  = bc ((XI - K + 0));
-        AY_1 = bc ((YI - L - 1));
-        AX_1 = bc ((XI - K - 1));
-        AY_2 = bc ((YI - L - 2));
-        AX_2 = bc ((XI - K - 2));
+        AY1  = bc (YI - L + 1);
+        AX1  = bc (XI - K + 1);
+        AY0  = bc (YI - L + 0);
+        AX0  = bc (XI - K + 0);
+        AY_1 = bc (YI - L - 1);
+        AX_1 = bc (XI - K - 1);
+        AY_2 = bc (YI - L - 2);
+        AX_2 = bc (XI - K - 2);
 
         ## Perform interpolation
         sz = size (Z);
-        ZI = AY_2 .* AX_2 .* Z (sym_sub2ind (sz, L+2, K+2)) ...
-           + AY_2 .* AX_1 .* Z (sym_sub2ind (sz, L+2, K+1)) ...
-           + AY_2 .* AX0  .* Z (sym_sub2ind (sz, L+2, K))   ...
-           + AY_2 .* AX1  .* Z (sym_sub2ind (sz, L+2, K-1)) ...
-           + AY_1 .* AX_2 .* Z (sym_sub2ind (sz, L+1, K+2)) ...
-           + AY_1 .* AX_1 .* Z (sym_sub2ind (sz, L+1, K+1)) ...
-           + AY_1 .* AX0  .* Z (sym_sub2ind (sz, L+1, K))   ...
-           + AY_1 .* AX1  .* Z (sym_sub2ind (sz, L+1, K-1)) ...
-           + AY0  .* AX_2 .* Z (sym_sub2ind (sz, L,   K+2)) ...
-           + AY0  .* AX_1 .* Z (sym_sub2ind (sz, L,   K+1)) ...
-           + AY0  .* AX0  .* Z (sym_sub2ind (sz, L,   K))   ...
-           + AY0  .* AX1  .* Z (sym_sub2ind (sz, L,   K-1)) ...
-           + AY1  .* AX_2 .* Z (sym_sub2ind (sz, L-1, K+2)) ...
-           + AY1  .* AX_1 .* Z (sym_sub2ind (sz, L-1, K+1)) ...
-           + AY1  .* AX0  .* Z (sym_sub2ind (sz, L-1, K))   ...
-           + AY1  .* AX1  .* Z (sym_sub2ind (sz, L-1, K-1));
-        ZI (!inside) = extrapval;
+        ZI = AY_2 .* AX_2 .* Z(sym_sub2ind (sz, L+2, K+2)) ...
+           + AY_2 .* AX_1 .* Z(sym_sub2ind (sz, L+2, K+1)) ...
+           + AY_2 .* AX0  .* Z(sym_sub2ind (sz, L+2, K))   ...
+           + AY_2 .* AX1  .* Z(sym_sub2ind (sz, L+2, K-1)) ...
+           + AY_1 .* AX_2 .* Z(sym_sub2ind (sz, L+1, K+2)) ...
+           + AY_1 .* AX_1 .* Z(sym_sub2ind (sz, L+1, K+1)) ...
+           + AY_1 .* AX0  .* Z(sym_sub2ind (sz, L+1, K))   ...
+           + AY_1 .* AX1  .* Z(sym_sub2ind (sz, L+1, K-1)) ...
+           + AY0  .* AX_2 .* Z(sym_sub2ind (sz, L,   K+2)) ...
+           + AY0  .* AX_1 .* Z(sym_sub2ind (sz, L,   K+1)) ...
+           + AY0  .* AX0  .* Z(sym_sub2ind (sz, L,   K))   ...
+           + AY0  .* AX1  .* Z(sym_sub2ind (sz, L,   K-1)) ...
+           + AY1  .* AX_2 .* Z(sym_sub2ind (sz, L-1, K+2)) ...
+           + AY1  .* AX_1 .* Z(sym_sub2ind (sz, L-1, K+1)) ...
+           + AY1  .* AX0  .* Z(sym_sub2ind (sz, L-1, K))   ...
+           + AY1  .* AX1  .* Z(sym_sub2ind (sz, L-1, K-1));
+        ZI (!inside) = extrap;
 
       else
         error ("interp2: input data must have 'meshgrid' format");
       endif
+    #}
 
-    elseif (strcmp (method, "spline"))
+    if (strcmp (method, "spline"))
       if (isgriddata (XI) && isgriddata (YI'))
-        ZI = __splinen__ ({Y(:,1).', X(1,:)}, Z, {YI(:,1), XI(1,:)}, extrapval,
-                        "spline");
+        ZI = __splinen__ ({Y(:,1).', X(1,:)}, Z, {YI(:,1), XI(1,:)}, extrap,
+                          "spline");
       else
-        error ("interp2: input data must have 'meshgrid' format");
+        error ("interp2: XI, YI must have uniform spacing ('meshgrid' format)");
       endif
     else
-      error ("interp2: interpolation METHOD not recognized");
+      error ("interp2: unrecognized interpolation method '%s'", method);
     endif
 
   endif
+
 endfunction
 
 function b = isgriddata (X)
   d1 = diff (X, 1, 1);
-  b = all (d1 (:) == 0);
+  b = all (d1(:) == 0);
 endfunction
 
 ## Compute the bicubic interpolation coefficients
@@ -437,16 +450,16 @@
 
 ## This version of sub2ind behaves as if the data was symmetrically padded
 function ind = sym_sub2ind (sz, Y, X)
-  Y (Y < 1) = 1 - Y (Y < 1);
-  while (any (Y (:) > 2 * sz (1)))
-    Y (Y > 2 * sz (1)) = round (Y (Y > 2 * sz (1)) / 2);
+  Y(Y < 1) = 1 - Y(Y < 1);
+  while (any (Y(:) > 2*sz(1)))
+    Y(Y > 2*sz(1)) = round (Y(Y > 2*sz(1)) / 2);
   endwhile
-  Y (Y > sz (1)) = 1 + 2 * sz (1) - Y (Y > sz (1));
-  X (X < 1) = 1 - X (X < 1);
-  while (any (X (:) > 2 * sz (2)))
-    X (X > 2 * sz (2)) = round (X (X > 2 * sz (2)) / 2);
+  Y(Y > sz(1)) = 1 + 2*sz(1) - Y(Y > sz(1));
+  X(X < 1) = 1 - X(X < 1);
+  while (any (X(:) > 2*sz(2)))
+    X(X > 2 * sz(2)) = round (X(X > 2*sz(2)) / 2);
   endwhile
-  X (X > sz (2)) = 1 + 2 * sz (2) - X (X > sz (2));
+  X(X > sz(2)) = 1 + 2*sz(2) - X(X > sz(2));
   ind = sub2ind (sz, Y, X);
 endfunction
 
@@ -460,7 +473,7 @@
 %! yi = linspace (min (y), max (y), 26)';
 %! mesh (xi,yi,interp2 (x,y,A,xi,yi, "linear"));
 %! [x,y] = meshgrid (x,y);
-%! hold on; plot3 (x(:),y(:),A(:),"b*"); hold off;
+%! hold on; plot3 (x,y,A,"b*"); hold off;
 
 %!demo
 %! clf;
@@ -471,7 +484,7 @@
 %! yi = linspace (min (y), max (y), 41)';
 %! mesh (xi,yi,interp2 (x,y,A,xi,yi, "linear"));
 %! [x,y] = meshgrid (x,y);
-%! hold on; plot3 (x(:),y(:),A(:),"b*"); hold off;
+%! hold on; plot3 (x,y,A,"b*"); hold off;
 
 %!demo
 %! clf;
@@ -482,7 +495,7 @@
 %! yi = linspace (min (y), max (y), 26)';
 %! mesh (xi,yi,interp2 (x,y,A,xi,yi, "nearest"));
 %! [x,y] = meshgrid (x,y);
-%! hold on; plot3 (x(:),y(:),A(:),"b*"); hold off;
+%! hold on; plot3 (x,y,A,"b*"); hold off;
 
 %!demo
 %! clf;
@@ -493,9 +506,10 @@
 %! yi = linspace (min (y), max (y), 41)';
 %! mesh (xi,yi,interp2 (x,y,A,xi,yi, "nearest"));
 %! [x,y] = meshgrid (x,y);
-%! hold on; plot3 (x(:),y(:),A(:),"b*"); hold off;
+%! hold on; plot3 (x,y,A,"b*"); hold off;
 
-%!demo
+## 'pchip' commented out since it is the same as 'cubic'
+%!#demo
 %! clf;
 %! colormap ("default");
 %! A = [13,-1,12;5,4,3;1,6,2];
@@ -504,9 +518,10 @@
 %! yi = linspace (min (y), max (y), 26)';
 %! mesh (xi,yi,interp2 (x,y,A,xi,yi, "pchip"));
 %! [x,y] = meshgrid (x,y);
-%! hold on; plot3 (x(:),y(:),A(:),"b*"); hold off;
+%! hold on; plot3 (x,y,A,"b*"); hold off;
 
-%!demo
+## 'pchip' commented out since it is the same as 'cubic'
+%!#demo
 %! clf;
 %! colormap ("default");
 %! [x,y,A] = peaks (10);
@@ -515,7 +530,7 @@
 %! yi = linspace (min (y), max (y), 41)';
 %! mesh (xi,yi,interp2 (x,y,A,xi,yi, "pchip"));
 %! [x,y] = meshgrid (x,y);
-%! hold on; plot3 (x(:),y(:),A(:),"b*"); hold off;
+%! hold on; plot3 (x,y,A,"b*"); hold off;
 
 %!demo
 %! clf;
@@ -526,7 +541,7 @@
 %! yi = linspace (min (y), max (y), 26)';
 %! mesh (xi,yi,interp2 (x,y,A,xi,yi, "cubic"));
 %! [x,y] = meshgrid (x,y);
-%! hold on; plot3 (x(:),y(:),A(:),"b*"); hold off;
+%! hold on; plot3 (x,y,A,"b*"); hold off;
 
 %!demo
 %! clf;
@@ -537,7 +552,7 @@
 %! yi = linspace (min (y), max (y), 41)';
 %! mesh (xi,yi,interp2 (x,y,A,xi,yi, "cubic"));
 %! [x,y] = meshgrid (x,y);
-%! hold on; plot3 (x(:),y(:),A(:),"b*"); hold off;
+%! hold on; plot3 (x,y,A,"b*"); hold off;
 
 %!demo
 %! clf;
@@ -548,7 +563,7 @@
 %! yi = linspace (min (y), max (y), 26)';
 %! mesh (xi,yi,interp2 (x,y,A,xi,yi, "spline"));
 %! [x,y] = meshgrid (x,y);
-%! hold on; plot3 (x(:),y(:),A(:),"b*"); hold off;
+%! hold on; plot3 (x,y,A,"b*"); hold off;
 
 %!demo
 %! clf;
@@ -559,61 +574,63 @@
 %! yi = linspace (min (y), max (y), 41)';
 %! mesh (xi,yi,interp2 (x,y,A,xi,yi, "spline"));
 %! [x,y] = meshgrid (x,y);
-%! hold on; plot3 (x(:),y(:),A(:),"b*"); hold off;
+%! hold on; plot3 (x,y,A,"b*"); hold off;
 
-%!test % simple test
+%!test  # simple test
 %! x = [1,2,3];
 %! y = [4,5,6,7];
 %! [X, Y] = meshgrid (x, y);
-%! Orig = X.^2 + Y.^3;
+%! orig = X.^2 + Y.^3;
 %! xi = [1.2,2, 1.5];
 %! yi = [6.2, 4.0, 5.0]';
 %!
-%! Expected = ...
+%! expected = ...
 %!   [243,   245.4,  243.9;
 %!     65.6,  68,     66.5;
 %!    126.6, 129,    127.5];
-%! Result = interp2 (x,y,Orig, xi, yi);
+%! result = interp2 (x,y,orig, xi, yi);
 %!
-%! assert (Result, Expected, 1000*eps);
+%! assert (result, expected, 1000*eps);
 
-%!test % 2^n form
+%!test  # 2^n refinement form
 %! x = [1,2,3];
 %! y = [4,5,6,7];
 %! [X, Y] = meshgrid (x, y);
-%! Orig = X.^2 + Y.^3;
+%! orig = X.^2 + Y.^3;
 %! xi = [1:0.25:3];  yi = [4:0.25:7]';
-%! Expected = interp2 (x,y,Orig, xi, yi);
-%! Result = interp2 (Orig, 2);
+%! expected = interp2 (x,y,orig, xi, yi);
+%! result = interp2 (orig, 2);
 %!
-%! assert (Result, Expected, 10*eps);
+%! assert (result, expected, 10*eps);
 
-%!test % matrix slice
+%!test  # matrix slice
 %! A = eye (4);
 %! assert (interp2 (A,[1:4],[1:4]), [1,1,1,1]);
 
-%!test % non-gridded XI,YI
+%!test  # non-gridded XI,YI
 %! A = eye (4);
 %! assert (interp2 (A,[1,2;3,4],[1,3;2,4]), [1,0;0,1]);
 
-%!test % for values outside of boundaries
+%!test  # for values outside of boundaries
 %! x = [1,2,3];
 %! y = [4,5,6,7];
 %! [X, Y] = meshgrid (x,y);
-%! Orig = X.^2 + Y.^3;
+%! orig = X.^2 + Y.^3;
 %! xi = [0,4];
 %! yi = [3,8]';
-%! assert (interp2 (x,y,Orig, xi, yi), [NA,NA;NA,NA]);
-%! assert (interp2 (x,y,Orig, xi, yi,"linear", 0), [0,0;0,0]);
+%! assert (interp2 (x,y,orig, xi, yi), [NA,NA;NA,NA]);
+%! assert (interp2 (x,y,orig, xi, yi,"linear", 0), [0,0;0,0]);
+%! assert (interp2 (x,y,orig, xi, yi,"linear", "extrap"), [1,17;468,484]);
+%! assert (interp2 (x,y,orig, xi, yi,"nearest", "extrap"), orig([1,end],[1,end]));
 
-%!test % for values at boundaries
-%! A=[1,2;3,4];
-%! x=[0,1];
-%! y=[2,3]';
+%!test  # for values at boundaries
+%! A = [1,2;3,4];
+%! x = [0,1];
+%! y = [2,3]';
 %! assert (interp2 (x,y,A,x,y,"linear"), A);
 %! assert (interp2 (x,y,A,x,y,"nearest"), A);
 
-%!test % for Matlab-compatible rounding for 'nearest'
+%!test  # for Matlab-compatible rounding for 'nearest'
 %! X = meshgrid (1:4);
 %! assert (interp2 (X, 2.5, 2.5, "nearest"), 3);
 
@@ -621,6 +638,7 @@
 %! z = [1 3 5; 3 5 7; 5 7 9];
 %! zout = [1 2 3 4 5; 2 3 4 5 6; 3 4 5 6 7; 4 5 6 7 8; 5 6 7 8 9];
 %! tol = 2 * eps;
+%!
 %!assert (interp2 (z), zout, tol)
 %!assert (interp2 (z, "linear"), zout, tol)
 %!assert (interp2 (z, "pchip"), zout, tol)
@@ -636,5 +654,34 @@
 %!assert (interp2 (z, [2 3 1], [2 2 2], "spline"), [5 7 3], tol)
 
 %% Test input validation
+%!error <Z must be a matrix> interp2 ({1})
+%!error <N must be an integer .= 0> interp2 (1, ones (2))
+%!error <N must be an integer .= 0> interp2 (1, -1)
+%!error <N must be an integer .= 0> interp2 (1, 1.5)
+%!error <METHOD must be a string> interp2 (1, 1, 1, 1, 1, 2)
+%!error <EXTRAP must be a numeric scalar or "extrap"> interp2 (1, 1, 1, 1, 1, 'linear', {1})
+%!error <EXTRAP must be a numeric scalar or "extrap"> interp2 (1, 1, 1, 1, 1, 'linear', ones (2,2))
+%!error <EXTRAP must be a numeric scalar or "extrap"> interp2 (1, 1, 1, 1, 1, 'linear', "abc")
 %!warning <ignoring unsupported '\*' flag> interp2 (rand (3,3), 1, "*linear");
+%!error <X, Y must be numeric matrices> interp2 ({1}, 1, 1, 1, 1)
+%!error <X, Y must be numeric matrices> interp2 (1, {1}, 1, 1, 1)
+%!error <XI, YI must be numeric> interp2 (1, 1, 1, {1}, 1)
+%!error <XI, YI must be numeric> interp2 (1, 1, 1, 1, {1})
+%!error <X and Y must be matrices of equal size> interp2 (ones(2,2), 1, 1, 1, 1)
+%!error <X and Y must be matrices of equal size> interp2 (ones(2,2), ones(2,3), 1, 1, 1)
+%!error <X and Y size must match the dimensions of Z> interp2 (1:3, 1:3, ones (3,2), 1, 1)
+%!error <X and Y size must match the dimensions of Z> interp2 (1:2, 1:2, ones (3,2), 1, 1)
+%!error <XI and YI must be matrices of equal size> interp2 (1, 1, 1, ones(2,2), 1)
+%!error <XI and YI must be matrices of equal size> interp2 (1, 1, 1, 1, ones(2,2))
+%!error <pchip requires at least 2 points> interp2 (1, 1, 1, 1, 1, "pchip")
+%!error <cubic requires at least 2 points> interp2 (1, 1, 1, 1, 1, "cubic")
+%!error <X and Y size must match the dimensions of Z> interp2 (1:3, 1:3, ones (3,2), 1, 1, "spline")
+%!error <X and Y size must match the dimensions of Z> interp2 (1:2, 1:2, ones (3,2), 1, 1, "spline")
+%!error <X and Y must be matrices of equal size> interp2 (ones(2,2), 1, 1, 1, 1, "spline")
+%!error <X and Y must be matrices of equal size> interp2 (ones(2,2), ones(2,3), 1, 1, 1, "spline")
+%!error <XI and YI must be matrices of equal size> interp2 (1, 1, 1, ones(2,2), 1, "spline")
+%!error <XI and YI must be matrices of equal size> interp2 (1, 1, 1, 1, ones(2,2), "spline")
+%!error <XI, YI must have uniform spacing> interp2 (1, 1, 1, [1 2 4], [1 2 3], "spline")
+%!error <XI, YI must have uniform spacing> interp2 (1, 1, 1, [1 2 3], [1 2 4], "spline")
+%!error <unrecognized interpolation method 'foobar'> interp2 (1, 1, 1, 1, 1, "foobar")
 

--- a/scripts/general/interp3.m	Sun Mar 30 13:02:08 2014 -0700
+++ b/scripts/general/interp3.m	Sun Mar 30 14:18:43 2014 -0700
@@ -19,52 +19,72 @@
 ## -*- texinfo -*-
 ## @deftypefn  {Function File} {@var{vi} =} interp3 (@var{x}, @var{y}, @var{z}, @var{v}, @var{xi}, @var{yi}, @var{zi})
 ## @deftypefnx {Function File} {@var{vi} =} interp3 (@var{v}, @var{xi}, @var{yi}, @var{zi})
-## @deftypefnx {Function File} {@var{vi} =} interp3 (@var{v}, @var{m})
+## @deftypefnx {Function File} {@var{vi} =} interp3 (@var{v}, @var{n})
 ## @deftypefnx {Function File} {@var{vi} =} interp3 (@var{v})
 ## @deftypefnx {Function File} {@var{vi} =} interp3 (@dots{}, @var{method})
 ## @deftypefnx {Function File} {@var{vi} =} interp3 (@dots{}, @var{method}, @var{extrapval})
 ##
-## Perform 3-dimensional interpolation.  Each element of the 3-dimensional
-## array @var{v} represents a value at a location given by the parameters
-## @var{x}, @var{y}, and @var{z}.  The parameters @var{x}, @var{x}, and
-## @var{z} are either 3-dimensional arrays of the same size as the array
-## @var{v} in the @qcode{"meshgrid"} format or vectors.  The parameters
-## @var{xi}, etc. respect a similar format to @var{x}, etc., and they
-## represent the points at which the array @var{vi} is interpolated.
+## Three-dimensional interpolation.
+##
+## Interpolate reference data @var{x}, @var{y}, @var{z}, @var{v} to determine
+## @var{vi} at the coordinates @var{xi}, @var{yi}, @var{zi}.  The reference
+## data @var{x}, @var{y}, @var{z} can be matrices, as returned by
+## @code{meshgrid}, in which case the sizes of
+## @var{x}, @var{y}, @var{z}, and @var{v} must be equal.  If @var{x}, @var{y},
+## @var{z} are vectors describing a cubic grid then
+## @code{length (@var{x}) == columns (@var{v})},
+## @code{length (@var{y}) == rows (@var{v})},
+## and @code{length (@var{z}) == size (@var{v}, 3)}.  In either case the input
+## data must be strictly monotonic.
 ##
-## If @var{x}, @var{y}, @var{z} are omitted, they are assumed to be
-## @code{x = 1 : size (@var{v}, 2)}, @code{y = 1 : size (@var{v}, 1)} and
-## @code{z = 1 : size (@var{v}, 3)}.  If @var{m} is specified, then
-## the interpolation adds a point half way between each of the interpolation
-## points.  This process is performed @var{m} times.  If only @var{v} is
-## specified, then @var{m} is assumed to be @code{1}.
+## If called without @var{x}, @var{y}, @var{z}, and just a single reference
+## data matrix @var{v}, the 3-D region
+## @code{@var{x} = 1:columns (@var{v}), @var{y} = 1:rows (@var{v}),
+## @var{z} = 1:size (@var{v}, 3)} is assumed.
+## This saves memory if the grid is regular and the distance between points is
+## not important.
 ##
-## Method is one of:
+## If called with a single reference data matrix @var{v} and a refinement
+## value @var{n}, then perform interpolation over a 3-D grid where each original
+## interval has been recursively subdivided @var{n} times.  This results in
+## @code{2^@var{n}-1} additional points for every interval in the original
+## grid.  If @var{n} is omitted a value of 1 is used.  As an example, the
+## interval [0,1] with @code{@var{n}==2} results in a refined interval with
+## points at [0, 1/4, 1/2, 3/4, 1].
+##
+## The interpolation @var{method} is one of:
 ##
 ## @table @asis
 ## @item @qcode{"nearest"}
 ## Return the nearest neighbor.
 ##
-## @item @qcode{"linear"}
+## @item @qcode{"linear"} (default)
 ## Linear interpolation from nearest neighbors.
 ##
+## @item @qcode{"pchip"}
+## Piecewise cubic Hermite interpolating polynomial---shape-preserving
+## interpolation with smooth first derivative (not implemented yet).
+##
 ## @item @qcode{"cubic"}
-## Cubic interpolation from four nearest neighbors (not implemented yet).
+## Cubic interpolation (same as @qcode{"pchip"} [not implemented yet]).
 ##
 ## @item @qcode{"spline"}
 ## Cubic spline interpolation---smooth first and second derivatives
 ## throughout the curve.
 ## @end table
 ##
-## The default method is @qcode{"linear"}.
-##
-## If @var{extrap} is the string @qcode{"extrap"}, then extrapolate values
-## beyond the endpoints.  If @var{extrap} is a number, replace values beyond
-## the endpoints with that number.  If @var{extrap} is missing, assume NA.
-## @seealso{interp1, interp2, spline, meshgrid}
+## If @var{extrapval} is a number, then replace values beyond the endpoints
+## with that number.  When unspecified, @var{extrapval} defaults to @code{NA}.
+## Note that if @var{extrapval} is used, @var{method} must be specified as well.
+## @seealso{interp1, interp2, interpn, meshgrid}
 ## @end deftypefn
 
+## FIXME: Need to validate N argument (maybe change interpn).
+## FIXME: Need to add support for 'pchip' method (maybe change interpn).
+## FIXME: Need to add support for "extrap" string value (maybe change interpn).
+
 function vi = interp3 (varargin)
+
   method = "linear";
   extrapval = NA;
   nargs = nargin;
@@ -75,14 +95,15 @@
 
   if (ischar (varargin{end}))
     method = varargin{end};
-    nargs = nargs - 1;
-  elseif (nargs > 1 && ischar (varargin{end - 1}))
+    nargs--;
+  elseif (nargs > 1 && ischar (varargin{end-1}))
+    ## FIXME: No support for "extrap" string
     if (! isnumeric (varargin{end}) || ! isscalar (varargin{end}))
-      error ("interp3: extrapal is expected to be a numeric scalar");
+      error ("interp3: EXTRAPVAL must be a numeric scalar");
     endif
     extrapval = varargin{end};
     method = varargin{end-1};
-    nargs = nargs - 2;
+    nargs -= 2;
   endif
 
   if (method(1) == "*")
@@ -90,102 +111,127 @@
     method(1) = [];
   endif
 
-  if (nargs < 3 || (nargs == 4 && ! isvector (varargin{1})
-                    && nargs == (ndims (varargin{1}) + 1)))
+  if (nargs < 3)
+    ## Calling form interp3 (v) OR interp3 (v, n)
+    v = varargin{1};
+    if (ndims (v) != 3)
+      error ("interp3: V must be a 3-D array of values");
+    endif
+    n = varargin(2:nargs);
+    v = permute (v, [2, 1, 3]);
+    vi = ipermute (interpn (v, n{:}, method, extrapval), [2, 1, 3]);
+
+  elseif (nargs == 4 && ! isvector (varargin{1}))
+    ## Calling form interp3 (v, xi, yi, zi)
     v = varargin{1};
     if (ndims (v) != 3)
-      error ("interp3: expect 3-dimensional array of values");
+      error ("interp3: V must be a 3-D array of values");
     endif
-    x = varargin (2:nargs);
-    if (any (! cellfun (@isvector, x)))
-      for i = 2 : 3
-        if (! size_equal (x{1}, x{i}))
-          error ("interp3: dimensional mismatch");
-        endif
-        x{i} = permute (x{i}, [2, 1, 3]);
+    xi = varargin(2:4);
+    if (any (! cellfun (@isvector, xi)))
+      ## Meshgridded values rather than vectors
+      if (! size_equal (xi{:}))
+        error ("interp3: XI, YI, and ZI dimensions must be equal");
+      endif
+      for i = 1 : 3
+        xi{i} = permute (xi{i}, [2, 1, 3]);
       endfor
-      x{1} = permute (x{1}, [2, 1, 3]);
     endif
     v = permute (v, [2, 1, 3]);
-    vi = ipermute (interpn (v, x{:}, method, extrapval), [2, 1, 3]);
-  elseif (nargs == 7 && nargs == (2 * ndims (varargin{ceil (nargs / 2)})) + 1)
+    vi = ipermute (interpn (v, xi{:}, method, extrapval), [2, 1, 3]);
+
+  elseif (nargs == 7)
+    ## Calling form interp3 (x, y, z, v, xi, yi, zi)
     v = varargin{4};
     if (ndims (v) != 3)
-      error ("interp3: expect 3-dimensional array of values");
+      error ("interp3: V must be a 3-D array of values");
     endif
-    x = varargin (1:3);
+    x = varargin(1:3);
     if (any (! cellfun (@isvector, x)))
-      for i = 2 : 3
-        if (! size_equal (x{1}, x{i}) || ! size_equal (x{i}, v))
-          error ("interp3: dimensional mismatch");
-        endif
+      ## Meshgridded values rather than vectors
+      if (! size_equal (x{:}, v))
+        error ("interp3: X, Y, Z, and V dimensions must be equal");
+      endif
+      for i = 1 : 3
         x{i} = permute (x{i}, [2, 1, 3]);
       endfor
-      x{1} = permute (x{1}, [2, 1, 3]);
     endif
-    y = varargin (5:7);
-    if (any (! cellfun (@isvector, y)))
-      for i = 2 : 3
-        if (! size_equal (y{1}, y{i}))
-          error ("interp3: dimensional mismatch");
-        endif
-        y{i} = permute (y{i}, [2, 1, 3]);
+    xi = varargin(5:7);
+    if (any (! cellfun (@isvector, xi)))
+      ## Meshgridded values rather than vectors
+      if (! size_equal (xi{:}))
+        error ("interp3: XI, YI, and ZI dimensions must be equal");
+      endif
+      for i = 1 : 3
+        xi{i} = permute (xi{i}, [2, 1, 3]);
       endfor
-      y{1} = permute (y{1}, [2, 1, 3]);
     endif
     v = permute (v, [2, 1, 3]);
-    vi = ipermute (interpn (x{:}, v, y{:}, method, extrapval), [2, 1, 3]);
+    vi = ipermute (interpn (x{:}, v, xi{:}, method, extrapval), [2, 1, 3]);
+
   else
     error ("interp3: wrong number or incorrectly formatted input arguments");
   endif
+
 endfunction
 
 
-%!test
-%! x = y = z = -1:1; y = y + 2;
+%% FIXME: Need some demo blocks here to show off the function like interp2.m.
+
+%!test  # basic test
+%! x = y = z = -1:1;  y = y + 2;
 %! f = @(x,y,z) x.^2 - y - z.^2;
 %! [xx, yy, zz] = meshgrid (x, y, z);
 %! v = f (xx,yy,zz);
-%! xi = yi = zi = -1:0.5:1; yi = yi + 2.1;
+%! xi = yi = zi = -1:0.5:1;  yi = yi + 2.1;
 %! [xxi, yyi, zzi] = meshgrid (xi, yi, zi);
 %! vi = interp3 (x, y, z, v, xxi, yyi, zzi);
 %! [xxi, yyi, zzi] = ndgrid (yi, xi, zi);
 %! vi2 = interpn (y, x, z, v, xxi, yyi, zzi);
-%! tol = 10 * eps;
-%! assert (vi, vi2, tol);
+%! assert (vi, vi2, 10*eps);
 
-%!test
-%! x=z=1:2; y=1:3;xi=zi=.6:1.6; yi=1; v=ones([3,2,2]);  v(:,2,1)=[7 ;5;4];  v(:,1,2)=[2 ;3;5];
+%!test  # meshgridded xi, yi, zi
+%! x = z = 1:2;  y = 1:3;
+%! v = ones ([3,2,2]);  v(:,2,1) = [7;5;4];  v(:,1,2) = [2;3;5];
+%! xi = zi = .6:1.6;  yi = 1; 
 %! [xxi3, yyi3, zzi3] = meshgrid (xi, yi, zi);
 %! [xxi, yyi, zzi] = ndgrid (yi, xi, zi);
 %! vi = interp3 (x, y, z, v, xxi3, yyi3, zzi3, "nearest");
-%! vi2 = interpn (y, x, z, v, xxi, yyi, zzi,"nearest");
+%! vi2 = interpn (y, x, z, v, xxi, yyi, zzi, "nearest");
 %! assert (vi, vi2);
 
-%!test
-%! x=z=1:2; y=1:3;xi=zi=.6:1.6; yi=1; v=ones([3,2,2]);  v(:,2,1)=[7 ;5;4];  v(:,1,2)=[2 ;3;5];
-%! vi = interp3 (x, y, z, v, xi+1, yi, zi, "nearest",3);
-%! vi2 = interpn (y, x, z, v, yi, xi+1, zi,"nearest", 3);
-%! assert (vi, vi2);
-
-%!test
-%! x=z=1:2; y=1:3;xi=zi=.6:1.6; yi=1; v=ones([3,2,2]);  v(:,2,1)=[7 ;5;4];  v(:,1,2)=[2 ;3;5];
+%!test  # vector xi, yi, zi
+%! x = z = 1:2;  y = 1:3;
+%! v = ones ([3,2,2]);  v(:,2,1) = [7;5;4];  v(:,1,2) = [2;3;5];
+%! xi = zi = .6:1.6;  yi = 1; 
 %! vi = interp3 (x, y, z, v, xi, yi, zi, "nearest");
 %! vi2 = interpn (y, x, z, v, yi, xi, zi,"nearest");
 %! assert (vi, vi2);
 
-%!test
-%! x=z=1:2; y=1:3;xi=zi=.6:1.6; yi=1; v=ones([3,2,2]);  v(:,2,1)=[7 ;5;4];  v(:,1,2)=[2 ;3;5];
-%! vi = interp3 (v, xi, yi, zi, "nearest",3);
-%! vi2 = interpn (v, yi, xi, zi,"nearest", 3);
+%!test  # vector xi+1 with extrap value
+%! x = z = 1:2;  y = 1:3;
+%! v = ones ([3,2,2]);  v(:,2,1) = [7;5;4];  v(:,1,2) = [2;3;5];
+%! xi = zi = .6:1.6;  yi = 1; 
+%! vi = interp3 (x, y, z, v, xi+1, yi, zi, "nearest", 3);
+%! vi2 = interpn (y, x, z, v, yi, xi+1, zi, "nearest", 3);
 %! assert (vi, vi2);
 
-%!test
-%! xi=zi=.6:1.6; yi=1; v=ones([3,2,2]);  v(:,2,1)=[7 ;5;4];  v(:,1,2)=[2 ;3;5];
+%!test  # input value matrix--no x,y,z
+%! x = z = 1:2;  y = 1:3;
+%! v = ones ([3,2,2]);  v(:,2,1) = [7;5;4];  v(:,1,2) = [2;3;5];
+%! xi = zi = .6:1.6;  yi = 1; 
 %! vi = interp3 (v, xi, yi, zi, "nearest");
 %! vi2 = interpn (v, yi, xi, zi,"nearest");
 %! assert (vi, vi2);
 
+%!test  # input value matrix--no x,y,z, with extrap value
+%! x = z = 1:2;  y = 1:3;
+%! v = ones ([3,2,2]);  v(:,2,1) = [7;5;4];  v(:,1,2) = [2;3;5];
+%! xi = zi = .6:1.6;  yi = 1; 
+%! vi = interp3 (v, xi, yi, zi, "nearest", 3);
+%! vi2 = interpn (v, yi, xi, zi, "nearest", 3);
+%! assert (vi, vi2);
+
 %!shared z, zout, tol
 %! z = zeros (3, 3, 3);
 %! zout = zeros (5, 5, 5);
@@ -200,9 +246,21 @@
 %!                  5 6 7 8 9] + (n-1);
 %! end
 %! tol = 10 * eps;
+%!
 %!assert (interp3 (z), zout, tol)
 %!assert (interp3 (z, "linear"), zout, tol)
 %!assert (interp3 (z, "spline"), zout, tol)
 
 %% Test input validation
+%!error interp3 ()
+%!error interp3 ({1})
+%!error <EXTRAPVAL must be a numeric scalar> interp3 (1,2,3,4,1,2,3,"linear", {1})
+%!error <EXTRAPVAL must be a numeric scalar> interp3 (1,2,3,4,1,2,3,"linear", ones (2,2))
 %!warning <ignoring unsupported '\*' flag> interp3 (rand (3,3,3), 1, "*linear");
+%!error <V must be a 3-D array> interp3 (rand (2,2))
+%!error <V must be a 3-D array> interp3 (rand (2,2), 1,1,1)
+%!error <XI, YI, and ZI dimensions must be equal> interp3 (rand (2,2,2), 1,1, ones (2,2))
+%!error <V must be a 3-D array> interp3 (1:2, 1:2, 1:2, rand (2,2), 1,1,1)
+%!error <X, Y, Z, and V dimensions must be equal> interp3 (ones(1,2,2), ones(2,2,2), ones(2,2,2), rand (2,2,2), 1,1,1)
+%!error <XI, YI, and ZI dimensions must be equal> interp3 (1:2, 1:2, 1:2, rand (2,2,2), 1,1, ones (2,2))
+