--- a/doc/ChangeLog
+++ b/doc/ChangeLog
@@ -1,3 +1,16 @@
+2007-06-12  David Bateman  <dbateman@free.fr>
+
+	* interpreter/interp.txi: Split into two section and document
+	interp3 and the differences in the treatement of the dimensions
+	between interpn and interp3.
+	* hashing.txi: Remove.
+	* system.txi: Move it here as a subsection. Include explanation
+	and example.
+	* interpreter/octave.texi: Add sections for the Interpolation
+	chapter. Remove references to Hashing chapter and hashing.texi,
+	and subsections for hashing to system utilities chapter.
+
+
 2007-06-12  2007-06-10  Søren Hauberg  <hauberg@gmail.com>
 
         * interpreter/diffeq.txi: Note that x-dot is the derivative of x.

--- a/doc/interpreter/hashing.txi
+++ b/doc/interpreter/hashing.txi
@@ -1,8 +0,0 @@
-@c Copyright (C) 2007 John W. Eaton
-@c This is part of the Octave manual.
-@c For copying conditions, see the file gpl.texi.
-
-@node Hashing Functions
-@chapter Hashing Functions
-
-@DOCSTRING(md5sum)

--- a/doc/interpreter/interp.txi
+++ b/doc/interpreter/interp.txi
@@ -5,18 +5,96 @@
 @node Interpolation
 @chapter Interpolation
 
-Octave provides the following functions for interpolation.
+@menu
+* One-dimensional Interpolation::
+* Multi-dimensional Interpolation::
+@end menu
+
+@node One-dimensional Interpolation
+@section One-dimensional Interpolation
 
 @DOCSTRING(interp1)
 
-@DOCSTRING(interp2)
+Fourier interpolation, is a resampling technique where a signal is
+converted to the frequency domain, padded with zeros and then
+reconverted to the time domain.
 
 @DOCSTRING(interpft)
 
-@DOCSTRING(interpn)
+There are two significant limitations on Fourier interpolation. Firstly,
+the function signal is assumed to be periodic, and so no periodic
+signals will be poorly represented at the edges. Secondly, both the
+signal and its interpolation are required to be sampled at equispaced
+points. An example of the use of @code{interpft} is
 
-@DOCSTRING(bicubic)
+@example
+@group
+t = 0 : 0.3 : pi; dt = t(2)-t(1);
+n = length (t); k = 100;
+ti = t(1) + [0 : k-1]*dt*n/k;
+y = sin (4*t + 0.3) .* cos (3*t - 0.1);
+yp = sin (4*ti + 0.3) .* cos (3*ti - 0.1);
+plot (ti, yp, 'g', ti, interp1(t, y, ti, 'spline'), 'b', ...
+      ti, interpft (y, k), 'c', t, y, 'r+');
+legend ('sin(4t+0.3)cos(3t-0.1','spline','interpft','data');
+@end group
+@end example
+
+which demonstrates the poor behavior of Fourier interpolation for non
+periodic functions.
+
+In additional the support function @code{spline} and @code{lookup} that
+underlie the @code{interp1} function can be called directly.
 
 @DOCSTRING(spline)
 
+The @code{lookup} is used by other interpolation function to identify
+the points of the original data that are closest to the current point
+of interest.
+
 @DOCSTRING(lookup)
+
+@node Multi-dimensional Interpolation
+@section Multi-dimensional Interpolation
+
+There are three multi-dimensional interpolation function in Octave, with
+similar capabilities.
+
+@DOCSTRING(interp2)
+
+@DOCSTRING(interp3)
+
+@DOCSTRING(interpn)
+
+A significant difference between @code{interpn} and the other two
+multidimensional interpolation function is the fashion in which the
+dimensions are treated. For @code{interp2} and @code{interp3}, the 'y'
+axis is considered to be the columns of the matrix, whereas the 'x'
+axis corresponds to the rows the the array. As Octave indexes arrays in
+column major order, the first dimension of any array is the columns, and
+so @code{interpn} effectively reverses the 'x' and 'y' dimensions. 
+Consider the example
+
+@example
+@group
+x = y = z = -1:1;
+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;
+[xxi, yyi, zzi] = meshgrid (xi, yi, zi);
+vi = interp3(x, y, z, v, xxi, yyi, zzi);
+[xxi, yyi, zzi] = ndgrid (xi, yi, zi);
+vi2 = interpn(x, y, z, v, xxi, yyi, zzi);
+@end group
+@end example
+
+@noindent
+where @code{vi} and @code{vi2} are identical. The reversal of the
+dimensions is treated in the @code{meshgrid} and @code{ndgrid} functions
+respectively.
+
+In additional the support function @code{bicubic} that underlies the
+cubic interpolation of @code{interp2} function can be called directly.
+
+@DOCSTRING(bicubic)

--- a/doc/interpreter/octave.texi
+++ b/doc/interpreter/octave.texi
@@ -158,7 +158,6 @@
 * Polynomial Manipulations::    
 * Interpolation::
 * Geometry::
-* Hashing Functions::
 * Control Theory::              
 * Signal Processing::           
 * Image Processing::            
@@ -437,7 +436,9 @@
 * Models::                      
 * Distributions::               
 
-Hashing Functions
+Interpolation
+* One-dimensional Interpolation::
+* Multi-dimensional Interpolation::
 
 Control Theory
 
@@ -488,6 +489,7 @@
 * Password Database Functions::  
 * Group Database Functions::    
 * System Information::          
+* Hashing Functions::
 
 Packages
 
@@ -587,7 +589,6 @@
 @include poly.texi
 @include interp.texi
 @include geometry.texi
-@include hashing.texi
 @include control.texi
 @include signal.texi
 @include image.texi

--- a/doc/interpreter/system.txi
+++ b/doc/interpreter/system.txi
@@ -23,6 +23,7 @@
 * Password Database Functions::  
 * Group Database Functions::    
 * System Information::          
+* Hashing Functions::
 @end menu
 
 @node Timing Utilities
@@ -398,3 +399,34 @@
 @DOCSTRING(octave_config_info)
 
 @DOCSTRING(getrusage)
+
+@node Hashing Functions
+@section Hashing Functions
+
+It is often necessary to find if two strings or files are
+identical. This might be done by comparing them character by character
+and looking for differences. However, this can be slow, and so comparing
+a hash of the string or file can be a rapid way of finding if the files
+differ.
+
+Another use of the hashing function is to check for file integrity. The
+user can check the hash of the file against a known value and find if
+the file they have is the same as the one that the original hash was
+produced with.
+
+Octave supplies the @code{md5sum} function to perfrom MD5 hashes on
+strings and files. An example of the use of @code{md5sum} function might
+be
+
+@example
+@group
+if exist (file, "file")
+  hash = md5sum (file);
+else
+  # Treat the variable "file" as a string
+  hash = md5sum (file, true);
+endif
+@end group
+@end example
+
+@DOCSTRING(md5sum)

--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,21 @@
+2007-06-12  David Bateman  <dbateman@free.fr>
+
+	* general/interp1.m: Change examples to use new graphics
+	interface.
+	* general/__splinen__.m: New support function for N-dimensional
+	spline interpolation.
+	* general/bicubic.m: Allow definition of extrapolation
+	value. Adapt tests to use new graphics interface
+	* general/interp2.m: Call __splinen__ for 2-D spline
+	interpolation. Make the lookup table code only be called for
+	linear and nearest methods.
+	* general/interpn.m: New function for N-dimensional, linear, nearest
+	and spline interpolation.
+	* general/interp3.m: New function for 3-dimensional, linear, nearest
+	and spline interpolation.
+	* polynomial/spline.m: Change examples to use new graphics
+	interface.
+	
 2007-06-12  Steve M. Robbins  <steve@sumost.ca>
 
 	* statistics/tests/wilcoxon_test.m: Error if N <= 25.

new file mode 100644
--- /dev/null
+++ b/scripts/general/__splinen__.m
@@ -0,0 +1,47 @@
+## Copyright (C) 2007 David Bateman
+##
+## 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 2, 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, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{yi} = } __splinen__ (@var{x}, @var{y}, @var{xi})
+## Internal support function for multi-dimensional splines.
+## @end deftypefn
+
+## FIXME: Allow arbitrary grids..
+
+function yi = __splinen__ (x, y, xi, extrapval, f)
+  if (nargin != 5)
+    error ("Incorrect number of arguments");
+  endif
+
+  if (!iscell (x) || length(x) < ndims(y) || any (! cellfun (@isvector, x)) ||
+      !iscell (xi) || length(xi) < ndims(y) || any (! cellfun (@isvector, xi)))
+    error ("%s: non gridded data or dimensions inconsistent", f);
+  endif
+  yi = y;
+  for i = length(x):-1:1
+    yi = spline (x{i}, yi, xi{i}).';
+  endfor
+
+  [xi{:}] = ndgrid (xi{:});
+  idx = zeros (size(xi{1}));
+  for i = 1 : length(x)
+    idx |= xi{i} < min (x{i}(:)) | xi{i} > max (x{i}(:));
+  endfor
+  yi(idx) = extrapval;
+endfunction

--- a/scripts/general/bicubic.m
+++ b/scripts/general/bicubic.m
@@ -18,32 +18,37 @@
 ## 02110-1301, USA.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {@var{zi}=} bicubic (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi})
+## @deftypefn {Function File} {@var{zi}=} bicubic (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}, @var{extrapval})
 ##
 ## Return a matrix @var{zi} corresponding to the bicubic
 ## interpolations at @var{xi} and @var{yi} of the data supplied
-## as @var{x}, @var{y} and @var{z}. 
+## as @var{x}, @var{y} and @var{z}. Points outside the grid are set
+## to @var{extrapval}
 ##
-## For further information please see bicubic.pdf available at
-## @url{http://wiki.woodpecker.org.cn/moin/Octave/Bicubic}
+## See @url{http://wiki.woodpecker.org.cn/moin/Octave/Bicubic}
+## for further information.
 ## @seealso{interp2}
 ## @end deftypefn
 
 ## Bicubic interpolation method.
 ## Author: Hoxide Ma <hoxide_dirac@yahoo.com.cn>
 
-function F = bicubic (X, Y, Z, XI, YI, spline_alpha)
+function F = bicubic (X, Y, Z, XI, YI, extrapval, spline_alpha)
 
-  if (nargin < 1 || nargin > 6)
+  if (nargin < 1 || nargin > 7)
     print_usage ();
   endif
 
-  if (nargin == 6 && prod (size (spline_alpha)) == 1)
+  if (nargin == 7 && isscalar(spline_alpha))
     a = spline_alpha
   else
     a = 0.5;
   endif
 
+  if (nargin < 6)
+    extrapval = NaN;
+  endif
+
   if (nargin <= 2)
     ## bicubic (Z) or bicubic (Z, 2)
     if (nargin == 1) 
@@ -177,12 +182,12 @@
 	 p(int,inds+2) .* cs2 + p(int,inds+3) .* cs3);
   endfor
 
-  ## set points outside the table to NaN
+  ## set points outside the table to extrapval
   if (! (isempty (xfirst_ind) && isempty (xlast_ind)))
-    F(:, [xfirst_ind, xlast_ind]) = NaN;
+    F(:, [xfirst_ind, xlast_ind]) = extrapval;
   endif
   if (! (isempty (yfirst_ind) && isempty (ylast_ind)))
-    F([yfirst_ind; ylast_ind], :) = NaN;
+    F([yfirst_ind; ylast_ind], :) = extrapval;
   endif
 
 endfunction
@@ -191,7 +196,7 @@
 %! A=[13,-1,12;5,4,3;1,6,2];
 %! x=[0,1,4]+10; y=[-10,-9,-8];
 %! xi=linspace(min(x),max(x),17);
-%! yi=linspace(min(y),max(y),26);
+%! yi=linspace(min(y),max(y),26)';
 %! mesh(xi,yi,bicubic(x,y,A,xi,yi));
 %! [x,y] = meshgrid(x,y);
 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;

--- a/scripts/general/interp1.m
+++ b/scripts/general/interp1.m
@@ -70,8 +70,9 @@
 ##    spl=interp1(xp,yp,xf,'spline');
 ##    cub=interp1(xp,yp,xf,'cubic');
 ##    near=interp1(xp,yp,xf,'nearest');
-##    plot(xf,yf,';original;',xf,lin,';linear;',xf,spl,';spline;',...
-##         xf,cub,';cubic;',xf,near,';nearest;',xp,yp,'*;;');
+##    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")
 ## @end group
 ## @end example
 ##
@@ -327,8 +328,8 @@
 %! spl=interp1(xp,yp,xf,"spline");
 %! cub=interp1(xp,yp,xf,"pchip");
 %! near=interp1(xp,yp,xf,"nearest");
-%! plot(xf,yf,";original;",xf,near,";nearest;",xf,lin,";linear;",...
-%!      xf,cub,";pchip;",xf,spl,";spline;",xp,yp,"*;;");
+%! 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")
 %! %--------------------------------------------------------
 %! % confirm that interpolated function matches the original
 
@@ -339,8 +340,8 @@
 %! spl=interp1(xp,yp,xf,"*spline");
 %! cub=interp1(xp,yp,xf,"*cubic");
 %! near=interp1(xp,yp,xf,"*nearest");
-%! plot(xf,yf,";*original;",xf,near,";*nearest;",xf,lin,";*linear;",...
-%!      xf,cub,";*cubic;",xf,spl,";*spline;",xp,yp,"*;;");
+%! 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")
 %! %--------------------------------------------------------
 %! % confirm that interpolated function matches the original
 

--- a/scripts/general/interp2.m
+++ b/scripts/general/interp2.m
@@ -62,7 +62,7 @@
 ## Cubic interpolation from four nearest neighbors.
 ## @item 'spline'
 ## Cubic spline interpolation--smooth first and second derivatives
-## throughout the curve (not implemented yet).
+## throughout the curve.
 ## @end table
 ##
 ## If a scalar value @var{extrapval} is defined as the final value, then
@@ -161,75 +161,102 @@
     error ("interp2 expected numeric XI, YI"); 
   endif
 
-  ## If X and Y vectors produce a grid from them
-  if (isvector (X) && isvector (Y))
-    [X, Y] = meshgrid (X, Y);
-  elseif (! size_equal (X, Y))
-    error ("X and Y must be matrices of same size");
-  endif
-  if (! size_equal (X, Z))
-    error ("X and Y size must match Z dimensions");
-  endif
+
+  if (strcmp (method, "linear") || strcmp (method, "nearest"))
+
+    ## If X and Y vectors produce a grid from them
+    if (isvector (X) && isvector (Y))
+      [X, Y] = meshgrid (X, Y);
+    elseif (! size_equal (X, Y))
+      error ("X and Y must be matrices of same size");
+    endif
+    if (! size_equal (X, Z))
+      error ("X and Y size must match Z dimensions");
+    endif
+
+    ## If Xi and Yi are vectors of different orientation build a grid
+    if ((rows (XI) == 1 && columns (YI) == 1)
+	|| (columns (XI) == 1 && rows (YI) == 1))
+      [XI, YI] = meshgrid (XI, YI);
+    elseif (! size_equal (XI, YI))
+      error ("XI and YI must be matrices of same size");
+    endif
 
-  ## If Xi and Yi are vectors of different orientation build a grid
-  if ((rows (XI) == 1 && columns (YI) == 1)
-      || (columns (XI) == 1 && rows (YI) == 1))
-    [XI, YI] = meshgrid (XI, YI);
-  elseif (! size_equal (XI, YI))
-    error ("XI and YI must be matrices of same size");
-  endif
+    shape = size (XI);
+    XI = reshape (XI, 1, prod (shape));
+    YI = reshape (YI, 1, prod (shape));
+
+    xidx = lookup (X(1, 2:end-1), XI) + 1;
+    yidx = lookup (Y(2:end-1, 1), YI) + 1;
 
-  shape = size (XI);
-  XI = reshape (XI, 1, prod (shape));
-  YI = reshape (YI, 1, prod (shape));
+    if (strcmp (method, "linear"))
+      ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy
+      ##
+      ## a-b
+      ## | |
+      ## c-d
+      a = Z(1:(zr - 1), 1:(zc - 1));
+      b = Z(1:(zr - 1), 2:zc) - a;
+      c = Z(2:zr, 1:(zc - 1)) - a;
+      d = Z(2:zr, 2:zc) - a - b - c;
 
-  xidx = lookup (X(1, 2:end-1), XI) + 1;
-  yidx = lookup (Y(2:end-1, 1), YI) + 1;
+      idx = sub2ind (size (a), yidx, xidx);
+
+      ## scale XI, YI values to a 1-spaced grid
+      Xsc = (XI - X(1, xidx)) ./ (X(1, xidx + 1) - X(1, xidx));
+      Ysc = (YI - Y(yidx, 1)') ./ (Y(yidx + 1, 1) - Y(yidx, 1))';
+
+      ## apply plane equation
+      ZI = a(idx) + b(idx).*Xsc + c(idx).*Ysc + d(idx).*Xsc.*Ysc;
 
-  if (strcmp (method, "linear"))
-    ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy
-    ##
-    ## a-b
-    ## | |
-    ## c-d
-    a = Z(1:(zr - 1), 1:(zc - 1));
-    b = Z(1:(zr - 1), 2:zc) - a;
-    c = Z(2:zr, 1:(zc - 1)) - a;
-    d = Z(2:zr, 2:zc) - a - b - c;
+    elseif (strcmp (method, "nearest"))
+      xtable = X(1, :);
+      ytable = Y(:, 1)';
+      ii = (XI - xtable(xidx) > xtable(xidx + 1) - XI);
+      jj = (YI - ytable(yidx) > ytable(yidx + 1) - YI);
+      idx = sub2ind (size (Z), yidx+jj, xidx+ii);
+      ZI = Z(idx);
+    endif
 
-    idx = sub2ind (size (a), yidx, xidx);
-
-    ## scale XI, YI values to a 1-spaced grid
-    Xsc = (XI - X(1, xidx)) ./ (X(1, xidx + 1) - X(1, xidx));
-    Ysc = (YI - Y(yidx, 1)') ./ (Y(yidx + 1, 1) - Y(yidx, 1))';
-
-    ## apply plane equation
-    ZI = a(idx) + b(idx).*Xsc + c(idx).*Ysc + d(idx).*Xsc.*Ysc;
+    ## set points outside the table to NaN
+    ZI (XI < X(1,1) | XI > X(1,end) | YI < Y(1,1) | YI > Y(end,1)) = extrapval;
+    ZI = reshape (ZI, shape);
+  else
 
-  elseif (strcmp (method, "nearest"))
-    xtable = X(1, :);
-    ytable = Y(:, 1)';
-    ii = (XI - xtable(xidx) > xtable(xidx + 1) - XI);
-    jj = (YI - ytable(yidx) > ytable(yidx + 1) - YI);
-    idx = sub2ind (size (Z), yidx+jj, xidx+ii);
-    ZI = Z(idx);
-
-  elseif (strcmp (method, "cubic"))
-    ## FIXME bicubic doesn't handle arbitrary XI, YI
-    ZI = bicubic (X, Y, Z, XI(1,:), YI(:,1));
+    ## If X and Y vectors produce a grid from them
+    if (isvector (X) && isvector (Y))
+      X = X(:).';
+      Y = Y(:);
+      if (!isequal ([length(X), length(Y)], size(Z)))
+	error ("X and Y size must match Z dimensions");
+      endif
+    elseif (!size_equal (X, Y))
+      error ("X and Y must be matrices of same size");
+      if (! size_equal (X, Z))
+	error ("X and Y size must match Z dimensions");
+      endif
+    endif
 
-  elseif (strcmp (method, "spline"))
-    ## FIXME Implement 2-D (or in fact ND) spline interpolation
-    error ("interp2, spline interpolation not yet implemented");
+    ## If Xi and Yi are vectors of different orientation build a grid
+    if ((rows (XI) == 1 && columns (YI) == 1)
+	|| (columns (XI) == 1 && rows (YI) == 1))
+      ## Do nothing
+    elseif (! size_equal (XI, YI))
+      error ("XI and YI must be matrices of same size");
+    endif
 
-  else
-    error ("interpolation method not recognized");
-  endif
+    ## FIXME bicubic/__splinen__ don't handle arbitrary XI, YI
+    if (strcmp (method, "cubic"))
+      ZI = bicubic (X, Y, Z, XI(1,:), YI(:,1), extrapval);
 
-  ## set points outside the table to NaN
-  ZI (XI < X(1,1) | XI > X(1,end) | YI < Y(1,1) | YI > Y(end,1)) = extrapval;
-  ZI = reshape (ZI, shape);
+    elseif (strcmp (method, "spline"))
+      ZI = __splinen__ ({Y(:,1).', X(1,:)}, Z, {YI(:,1), XI(1,:)}, extrapval, 
+			"spline");
+    else
+      error ("interpolation method not recognized");
+    endif
 
+  endif
 endfunction
 
 %!demo
@@ -250,15 +277,24 @@
 %! [x,y] = meshgrid(x,y); 
 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
 
-%!#demo
+%!demo
 %! A=[13,-1,12;5,4,3;1,6,2];
 %! x=[0,1,2]; y=[10,11,12];
 %! xi=linspace(min(x),max(x),17);
-%! yi=linspace(min(y),max(y),26);
+%! 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;
 
+%!demo
+%! A=[13,-1,12;5,4,3;1,6,2];
+%! x=[0,1,2]; y=[10,11,12];
+%! xi=linspace(min(x),max(x),17);
+%! 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;
+
 %!test % simple test
 %!  x = [1,2,3];
 %!  y = [4,5,6,7];

new file mode 100644
--- /dev/null
+++ b/scripts/general/interp3.m
@@ -0,0 +1,114 @@
+## Copyright (C) 2007 David Bateman
+##
+## 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 2, 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, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- 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})
+## @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 then 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 '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.
+##
+## If @var{x}, @var{y}, @var{z} are ommitted, 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 interplation 
+## points. This process is performed @var{m} times. If only @var{v} is 
+## specified, then @var{m} is assumed to be @code{1}.
+##
+## Method is one of:
+##
+## @table @asis
+## @item 'nearest'
+## Return the nearest neighbour.
+## @item 'linear'
+## Linear interpolation from nearest neighbours.
+## @item 'cubic'
+## Cubic interpolation from four nearest neighbours (not implemented yet).
+## @item 'spline'
+## Cubic spline interpolation--smooth first and second derivatives
+## throughout the curve.
+## @end table
+##
+## The default method is 'linear'.
+##
+## If @var{extrap} is the string '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 NaN.
+## @seealso{interp1, interp2, spline, meshgrid}
+## @end deftypefn
+
+function vi = interp3 (varargin)
+  method = "linear";
+  extrapval = NaN;
+  nargs = nargin;
+
+  if (nargin < 1)
+    print_usage ();
+  endif
+
+  if (ischar (varargin {end}))
+    method = varargin {end};
+    nargs = nargs - 1;
+  elseif (ischar (varargin {end - 1}))
+    if (! isnumeric (vargin {end}) || ! isscalar (vargin {end}))
+      error ("extrapal is expected to be a numeric scalar");
+    endif
+    method = varargin {end - 1};
+    nargs = nargs - 2;
+  endif
+
+  if (nargs < 3 || (nargs == 4 && ! isvector (varargin {1}) && 
+      nargs == (ndims (varargin {1}) + 1)))
+    v = varargin {1};
+    if (ndims (v) != 3)
+      error ("expect 3-dimensional array of values");
+    endif
+    vi = ipermute (interpn (permute(varargin, [1, 3, 2, 4]){:}), [2, 1, 3]);
+  elseif (nargs == 7 && nargs == (2 * ndims (varargin {ceil (nargs / 2)})) + 1)
+    v = varargin {4};
+    if (ndims (v) != 3)
+      error ("expect 3-dimensional array of values");
+    endif
+    vi = ipermute (interpn (permute(varargin, [2, 1, 3, 4, 6, 5, 7]){:}), 
+		   [2, 1, 3]);
+  else
+    error ("wrong number or incorrectly formatted input arguments");
+  endif
+endfunction
+
+%!test
+%! x = y = z = -1:1;
+%! 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;
+%! [xxi, yyi, zzi] = meshgrid (xi, yi, zi);
+%! vi = interp3(x, y, z, v, xxi, yyi, zzi);
+%! [xxi, yyi, zzi] = ndgrid (xi, yi, zi);
+%! vi2 = interpn(x, y, z, v, xxi, yyi, zzi);
+%! assert (vi, vi2);

--- a/scripts/general/interpft.m
+++ b/scripts/general/interpft.m
@@ -101,10 +101,9 @@
 %! ti = t(1) + [0 : k-1]*dt*n/k;
 %! y = sin (4*t + 0.3) .* cos (3*t - 0.1);
 %! yp = sin (4*ti + 0.3) .* cos (3*ti - 0.1);
-%! plot (ti, yp, 'g;sin(4t+0.3)cos(3t-0.1);', ...
-%!       ti, interp1 (t, y, ti, 'spline'), 'b;spline;', ...
-%!       ti, interpft (y ,k), 'c;interpft;', ...
-%!       t, y, 'r+;data;');
+%! plot (ti, yp, 'g', ti, interp1(t, y, ti, 'spline'), 'b', ...
+%!       ti, interpft (y, k), 'c', t, y, 'r+');
+%! legend ('sin(4t+0.3)cos(3t-0.1','spline','interpft','data');
 
 %!shared n,y
 %! x = [0:10]'; y = sin(x); n = length (x);

new file mode 100644
--- /dev/null
+++ b/scripts/general/interpn.m
@@ -0,0 +1,218 @@
+## Copyright (C) 2007 David Bateman
+##
+## 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 2, 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, write to the Free
+## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+## 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{vi} =} interpn (@var{x1}, @var{x2}, @dots{}, @var{v}, @var{y1}, @var{y2}, @dots{})
+## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v}, @var{y1}, @var{y2}, @dots{})
+## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v}, @var{m})
+## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v})
+## @deftypefnx {Function File} {@var{vi} =} interpn (@dots{}, @var{method})
+## @deftypefnx {Function File} {@var{vi} =} interpn (@dots{}, @var{method}, @var{extrapval})
+##
+## Perform @var{n}-dimensional interpolation, where @var{n} is at least two. 
+## Each element of then @var{n}-dimensional array @var{v} represents a value 
+## at a location given by the parameters @var{x1}, @var{x2}, @dots{}, @var{xn}. 
+## The parameters @var{x1}, @var{x2}, @dots{}, @var{xn} are either 
+## @var{n}-dimensional arrays of the same size as the array @var{v} in 
+## the 'ndgrid' format or vectors. The parameters @var{y1}, etc respect a 
+## similar format to @var{x1}, etc, and they represent the points at which
+## the array @var{vi} is interpolated.
+##
+## If @var{x1}, @dots{}, @var{xn} are ommitted, they are assumed to be 
+## @code{x1 = 1 : size (@var{v}, 1)}, etc. If @var{m} is specified, then
+## the interpolation adds a point half way between each of the interplation 
+## points. This process is performed @var{m} times. If only @var{v} is 
+## specified, then @var{m} is assumed to be @code{1}.
+##
+## Method is one of:
+##
+## @table @asis
+## @item 'nearest'
+## Return the nearest neighbour.
+## @item 'linear'
+## Linear interpolation from nearest neighbours.
+## @item 'cubic'
+## Cubic interpolation from four nearest neighbours (not implemented yet).
+## @item 'spline'
+## Cubic spline interpolation--smooth first and second derivatives
+## throughout the curve.
+## @end table
+##
+## The default method is 'linear'.
+##
+## If @var{extrap} is the string '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 NaN.
+## @seealso{interp1, interp2, spline, ndgrid}
+## @end deftypefn
+
+function vi = interpn (varargin)
+
+  method = "linear";
+  extrapval = NaN;
+  nargs = nargin;
+
+  if (nargin < 1)
+    print_usage ();
+  endif
+
+  if (ischar (varargin {end}))
+    method = varargin {end};
+    nargs = nargs - 1;
+  elseif (ischar (varargin {end - 1}))
+    if (! isnumeric (vargin {end}) || ! isscalar (vargin {end}))
+      error ("extrapal is expected to be a numeric scalar");
+    endif
+    method = varargin {end - 1};
+    nargs = nargs - 2;
+  endif
+
+  if (nargs < 3)
+    v = varargin {1};
+    m = 1;
+    if (nargs == 2)
+      m = varargin {2};
+      if (! isnumeric (m) || ! isscalar (m) || floor (m) != m)
+	error ("m is expected to be a integer scalar");
+      endif
+    endif
+    sz = size (v);
+    nd = ndims (v);
+    x = cell (1, nd);
+    y = cell (1, nd);
+    for i = 1 : nd;
+      x{i} = 1 : sz(i);
+      y{i} = 1 : (1 / (2 ^ m)) : sz(i);
+    endfor
+  elseif (! isvector (varargin {1}) && nargs == (ndims (varargin {1}) + 1))
+    v = varargin {1};
+    sz = size (v);
+    nd = ndims (v);
+    x = cell (1, nd);
+    y = varargin (2 : nargs);
+    for i = 1 : nd;
+      x{i} = 1 : sz(i);
+    endfor
+  elseif (rem (nargs, 2) == 1 && nargs ==  
+	  (2 * ndims (varargin {ceil (nargs / 2)})) + 1)
+    nv = ceil (nargs / 2);
+    v = varargin {nv};
+    sz = size (v);
+    nd = ndims (v);
+    x = varargin (1 : (nv - 1));
+    y = varargin ((nv + 1) : nargs);
+  else
+    error ("wrong number or incorrectly formatted input arguments");
+  endif
+
+  if (any (! cellfun (@isvector, x)))
+    for i = 2 : nd
+      if (! size_equal (x{1}, x{i}) || ! size_equal (x{i}, v))
+	error ("dimensional mismatch");
+      endif
+      idx (1 : nd) = {1};
+      idx (i) = ":";
+      x{i} = x{i}(idx{:});
+    endfor
+    idx (1 : nd) = {1};
+    idx (1) = ":";
+    x{1} = x{1}(idx{:});
+  endif
+
+  if (strcmp (method, "linear") || strcmp (method, "nearest"))
+    if (all (cellfun (@isvector, y)))
+      [y{:}] = ndgrid (y{:});
+    endif
+  elseif (any (! cellfun (@isvector, x)))
+    for i = 1 : nd
+      idx (1 : nd) = {1};
+      idx (i) = ":";
+      y{i} = y{i}(idx{:});
+    endfor
+  endif
+
+  method = tolower (method);
+  if (strcmp (method, "linear"))
+    vi = __lin_interpn__ (x{:}, v, y{:});
+    vi (vi == NaN) = extrapval;
+  elseif (strcmp (method, "nearest"))
+    yshape = size (y{1});
+    yidx = cell (1, nd);
+    for i = 1 : nd
+      y{i} = y{i}(:);
+      yidx{i} = lookup (x{i}(2:end-1), y{i}) + 1;
+    endfor
+    idx = cell (1,nd);
+    for i = 1 : nd
+      idx {i} = yidx{i} + (y{i} - x{i}(yidx{i}).' > ...
+			   x{i}(yidx{i} + 1).' - y{i});
+    endfor
+    vi = v (sub2ind (sz, idx{:}));
+    idx = zeros (prod(yshape),1);
+    for i = 1 : nd
+      idx |= y{i} < min (x{i}(:)) | y{i} > max (x{i}(:));
+    endfor
+    vi(idx) = extrapval;
+    vi = reshape (vi, yshape); 
+  elseif (strcmp (method, "spline")) 
+    vi = __splinen__ (x, v, y, extrapval, "interpn");
+  elseif (strcmp (method, "cubic")) 
+    error ("cubic interpolation not yet implemented");
+  else
+    error ("unrecognized interpolation method");
+  endif
+
+endfunction
+
+%!demo
+%! A=[13,-1,12;5,4,3;1,6,2];
+%! x=[0,1,4]; y=[10,11,12];
+%! xi=linspace(min(x),max(x),17);
+%! yi=linspace(min(y),max(y),26)';
+%! mesh(xi,yi,interpn(x,y,A.',xi,yi,"linear").');
+%! [x,y] = meshgrid(x,y); 
+%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
+
+%!demo
+%! A=[13,-1,12;5,4,3;1,6,2];
+%! x=[0,1,4]; y=[10,11,12];
+%! xi=linspace(min(x),max(x),17);
+%! yi=linspace(min(y),max(y),26)';
+%! mesh(xi,yi,interpn(x,y,A.',xi,yi,"nearest").');
+%! [x,y] = meshgrid(x,y); 
+%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
+
+%!#demo
+%! A=[13,-1,12;5,4,3;1,6,2];
+%! x=[0,1,2]; y=[10,11,12];
+%! xi=linspace(min(x),max(x),17);
+%! yi=linspace(min(y),max(y),26)';
+%! mesh(xi,yi,interpn(x,y,A.',xi,yi,"cubic").');
+%! [x,y] = meshgrid(x,y); 
+%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
+
+%!demo
+%! A=[13,-1,12;5,4,3;1,6,2];
+%! x=[0,1,2]; y=[10,11,12];
+%! xi=linspace(min(x),max(x),17);
+%! yi=linspace(min(y),max(y),26)';
+%! mesh(xi,yi,interpn(x,y,A.',xi,yi,"spline").');
+%! [x,y] = meshgrid(x,y); 
+%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
+

--- a/scripts/polynomial/spline.m
+++ b/scripts/polynomial/spline.m
@@ -212,9 +212,8 @@
 %! x = 0:10; y = sin(x);
 %! xspline = 0:0.1:10; yspline = spline(x,y,xspline);
 %! title("spline fit to points from sin(x)");
-%! plot(xspline,sin(xspline),";original;",...
-%!      xspline,yspline,"-;interpolation;",...
-%!      x,y,"+;interpolation points;");
+%! plot(xspline,sin(xspline),"r",xspline,yspline,"g-",x,y,"b+");
+%! legend("original","interpolation","interpolation points");
 %! %--------------------------------------------------------
 %! % confirm that interpolated function matches the original
 

--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,10 @@
+2007-06-12  David Bateman  <dbateman@free.fr>
+
+	* DLD-FUNCTIONS/interpn.cc: Remove it.
+	* DLD-FUNCTIONS/__lin_interpn__.cc: Move it. This is now a support
+	function of interpn.m.
+	* Makefile.in (DLD_XSRC): Remove interpn.cc and add __lin_interpn__.cc.
+
 2007-06-07  David Bateman  <dbateman@free.fr>
 
 	* ov-fcn-handles.cc (octave_fcn_handle::save_hdf5): More care that

new file mode 100644
--- /dev/null
+++ b/src/DLD-FUNCTIONS/__lin_interpn__.cc
@@ -0,0 +1,261 @@
+/*
+
+Copyright (C) 2007 Alexander Barth
+
+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 2, 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, write to the Free
+Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "dNDArray.h"
+
+#include "defun-dld.h"
+#include "error.h"
+#include "oct-obj.h"
+
+// equivalent to isvector.m
+
+bool
+isvector (const NDArray& array)
+{
+  const dim_vector dv = array.dims ();
+  return dv.length () == 2 && (dv(0) == 1 || dv(1) == 1);
+}
+
+// lookup a value in a sorted table (lookup.m)
+octave_idx_type
+lookup (const double *x, octave_idx_type n, double y)
+{
+  octave_idx_type j, j0, j1;
+
+  if (y > x[n-1] || y < x[0])
+    return -1;
+
+#ifdef EXHAUSTIF
+  for (j = 0; j < n - 1; j++)
+    {
+      if (x[j] <= y && y <= x[j+1])
+	return j;
+    }
+#else
+  j0 = 0;
+  j1 = n - 1;
+
+  while (true)
+    {
+      j = (j0+j1)/2;
+
+      if (y <= x[j+1])
+	{
+	  if (x[j] <= y)
+	    return j;
+
+	  j1 = j;
+	}
+
+      if (x[j] <= y)
+	j0 = j;
+    }
+
+#endif
+}
+
+// n-dimensional linear interpolation
+
+void
+lin_interpn (int n, const octave_idx_type *size, const octave_idx_type *scale,
+	     octave_idx_type Ni, double extrapval, const double **x,
+	     const double *v, const double **y, double *vi)
+{
+  bool out = false;
+  int bit;
+
+  OCTAVE_LOCAL_BUFFER (double, coef, 2*n);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, index, n);
+
+  // loop over all points
+  for (octave_idx_type m = 0; m < Ni; m++)
+    {
+      // loop over all dimensions
+      for (int i = 0; i < n; i++)
+	{
+          index[i] = lookup (x[i], size[i], y[i][m]);
+	  out = index[i] == -1;
+
+	  if (out)
+	    break;
+	  else
+            {
+	      octave_idx_type j = index[i];
+	      coef[2*i+1] = (y[i][m] - x[i][j])/(x[i][j+1] - x[i][j]);
+	      coef[2*i] = 1 - coef[2*i+1];
+	    }
+	}
+
+
+      if (out)
+	vi[m] = extrapval;
+      else
+	{
+	  vi[m] = 0;
+
+	  // loop over all corners of hypercube (1<<n = 2^n)
+	  for (int i = 0; i < (1 << n); i++)
+	    {
+	      double c = 1;
+	      octave_idx_type l = 0;
+
+	      // loop over all dimensions
+	      for (int j = 0; j < n; j++)
+		{
+		  // test if the jth bit in i is set
+		  bit = i >> j & 1;
+		  l += scale[j] * (index[j] + bit);
+		  c *= coef[2*j+bit];
+		}
+
+	      vi[m] += c * v[l];
+	    }
+	}
+    }
+}
+
+DEFUN_DLD (__lin_interpn__, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {@var{vi} =} __lin_interpn__ (@var{x1}, @var{x2}, @dots{}, @var{xn}, @var{v}, @var{y1}, @var{y2}, @dots{}, @var{yn})\n\
+Perform @var{n}-dimensional interpolation.  Each element of then\n\
+@var{n}-dimensional array @var{v} represents a value at a location\n\
+given by the parameters @var{x1}, @var{x2},...,@var{xn}. The parameters\n\
+@var{x1}, @var{x2}, @dots{}, @var{xn} are either @var{n}-dimensional\n\
+arrays of the same size as the array @var{v} in the \"ndgrid\" format\n\
+or vectors.  The parameters @var{y1}, @var{y2}, @dots{}, @var{yn} are\n\
+all @var{n}-dimensional arrays of the same size and represent the\n\
+points at which the array @var{vi} is interpolated.\n\
+\n\
+This function only performs linear interpolation.\n\
+@seealso{interp1, interp2, ndgrid}\n\
+@end deftypefn")
+{
+  octave_value retval;
+
+  int nargin = args.length ();
+
+  if (nargin < 2 ||  nargin % 2 == 0)
+    {
+      print_usage ();
+      return retval;
+    }
+
+  // dimension of the problem
+  int n = (nargin-1)/2;
+
+  OCTAVE_LOCAL_BUFFER (NDArray, X, n);
+  OCTAVE_LOCAL_BUFFER (NDArray, Y, n);
+
+  OCTAVE_LOCAL_BUFFER (const double *, x, n);
+  OCTAVE_LOCAL_BUFFER (const double *, y, n);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, scale, n);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, size, n);
+
+  const NDArray V = args(n).array_value ();
+  NDArray Vi = NDArray (args(n+1).dims ());
+
+  if (error_state)
+    {
+      print_usage ();
+      return retval;
+    }
+
+  const double *v = V.data ();
+  double *vi = Vi.fortran_vec ();
+  octave_idx_type Ni = Vi.numel ();
+
+  double extrapval = octave_NaN;
+
+  for (int i = 0; i < n; i++)
+    {
+      X[i] = args(i).array_value ();
+      Y[i] = args(n+i+1).array_value ();
+
+      if (error_state)
+	{
+	  print_usage ();
+	  return retval;
+	}
+
+      y[i] = Y[i].data ();
+      size[i] =  V.dims()(i);
+
+      if (Y[0].dims () != Y[i].dims ())
+	{
+	  error ("interpn: incompatible size of argument number %d", n+i+2);
+	  return retval;
+	}
+    }
+
+  // offset in memory of each dimension
+
+  scale[0] = 1;
+
+  for (int i = 1; i < n; i++)
+    scale[i] = scale[i-1] * size[i-1];
+
+  // tests if X[0] is a vector, if yes, assume that all elements of X are
+  // in the ndgrid format.
+
+  if (! isvector (X[0]))
+    {
+      for (int i = 0; i < n; i++)
+	{
+	  if (X[i].dims () != V.dims ())
+	    {
+	      error ("interpn: incompatible size of argument number %d", i+1);
+	      return retval;
+	    }
+	  else
+	    {
+              NDArray tmp = NDArray (dim_vector (size[i], 1));
+
+	      for (octave_idx_type j = 0; j < size[i]; j++)
+		tmp(j) =  X[i](scale[i]*j);
+
+              X[i] = tmp;
+	    }
+	}
+    }
+
+  for (int i = 0; i < n; i++)
+    {
+      if (! isvector (X[i]) && X[i].numel () != size[i])
+	{
+	  error ("interpn: incompatible size of argument number %d", i+1);
+	  return retval;
+	}
+      else
+	x[i] = X[i].data ();
+    }
+
+  lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi);
+
+  retval = Vi;
+
+  return retval;
+}

deleted file mode 100644
--- a/src/DLD-FUNCTIONS/interpn.cc
+++ /dev/null
@@ -1,261 +0,0 @@
-/*
-
-Copyright (C) 2007 Alexander Barth
-
-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 2, 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, write to the Free
-Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
-02110-1301, USA.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "dNDArray.h"
-
-#include "defun-dld.h"
-#include "error.h"
-#include "oct-obj.h"
-
-// equivalent to isvector.m
-
-bool
-isvector (const NDArray& array)
-{
-  const dim_vector dv = array.dims ();
-  return dv.length () == 2 && (dv(0) == 1 || dv(1) == 1);
-}
-
-// lookup a value in a sorted table (lookup.m)
-octave_idx_type
-lookup (const double *x, octave_idx_type n, double y)
-{
-  octave_idx_type j, j0, j1;
-
-  if (y > x[n-1] || y < x[0])
-    return -1;
-
-#ifdef EXHAUSTIF
-  for (j = 0; j < n - 1; j++)
-    {
-      if (x[j] <= y && y <= x[j+1])
-	return j;
-    }
-#else
-  j0 = 0;
-  j1 = n - 1;
-
-  while (true)
-    {
-      j = (j0+j1)/2;
-
-      if (y <= x[j+1])
-	{
-	  if (x[j] <= y)
-	    return j;
-
-	  j1 = j;
-	}
-
-      if (x[j] <= y)
-	j0 = j;
-    }
-
-#endif
-}
-
-// n-dimensional linear interpolation
-
-void
-lin_interpn (int n, const octave_idx_type *size, const octave_idx_type *scale,
-	     octave_idx_type Ni, double extrapval, const double **x,
-	     const double *v, const double **y, double *vi)
-{
-  bool out = false;
-  int bit;
-
-  OCTAVE_LOCAL_BUFFER (double, coef, 2*n);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, index, n);
-
-  // loop over all points
-  for (octave_idx_type m = 0; m < Ni; m++)
-    {
-      // loop over all dimensions
-      for (int i = 0; i < n; i++)
-	{
-          index[i] = lookup (x[i], size[i], y[i][m]);
-	  out = index[i] == -1;
-
-	  if (out)
-	    break;
-	  else
-            {
-	      octave_idx_type j = index[i];
-	      coef[2*i+1] = (y[i][m] - x[i][j])/(x[i][j+1] - x[i][j]);
-	      coef[2*i] = 1 - coef[2*i+1];
-	    }
-	}
-
-
-      if (out)
-	vi[m] = extrapval;
-      else
-	{
-	  vi[m] = 0;
-
-	  // loop over all corners of hypercube (1<<n = 2^n)
-	  for (int i = 0; i < (1 << n); i++)
-	    {
-	      double c = 1;
-	      octave_idx_type l = 0;
-
-	      // loop over all dimensions
-	      for (int j = 0; j < n; j++)
-		{
-		  // test if the jth bit in i is set
-		  bit = i >> j & 1;
-		  l += scale[j] * (index[j] + bit);
-		  c *= coef[2*j+bit];
-		}
-
-	      vi[m] += c * v[l];
-	    }
-	}
-    }
-}
-
-DEFUN_DLD (interpn, args, ,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {@var{vi} =} interpn (@var{x1}, @var{x2}, @dots{}, @var{xn}, @var{v}, @var{y1}, @var{y2}, @dots{}, @var{yn})\n\
-Perform @var{n}-dimensional interpolation.  Each element of then\n\
-@var{n}-dimensional array @var{v} represents a value at a location\n\
-given by the parameters @var{x1}, @var{x2},...,@var{xn}. The parameters\n\
-@var{x1}, @var{x2}, @dots{}, @var{xn} are either @var{n}-dimensional\n\
-arrays of the same size as the array @var{v} in the \"ndgrid\" format\n\
-or vectors.  The parameters @var{y1}, @var{y2}, @dots{}, @var{yn} are\n\
-all @var{n}-dimensional arrays of the same size and represent the\n\
-points at which the array @var{vi} is interpolated.\n\
-\n\
-This function only performs linear interpolation.\n\
-@seealso{interp1, interp2, ndgrid}\n\
-@end deftypefn")
-{
-  octave_value retval;
-
-  int nargin = args.length ();
-
-  if (nargin < 2 ||  nargin % 2 == 0)
-    {
-      print_usage ();
-      return retval;
-    }
-
-  // dimension of the problem
-  int n = (nargin-1)/2;
-
-  OCTAVE_LOCAL_BUFFER (NDArray, X, n);
-  OCTAVE_LOCAL_BUFFER (NDArray, Y, n);
-
-  OCTAVE_LOCAL_BUFFER (const double *, x, n);
-  OCTAVE_LOCAL_BUFFER (const double *, y, n);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, scale, n);
-  OCTAVE_LOCAL_BUFFER (octave_idx_type, size, n);
-
-  const NDArray V = args(n).array_value ();
-  NDArray Vi = NDArray (args(n+1).dims ());
-
-  if (error_state)
-    {
-      print_usage ();
-      return retval;
-    }
-
-  const double *v = V.data ();
-  double *vi = Vi.fortran_vec ();
-  octave_idx_type Ni = Vi.numel ();
-
-  double extrapval = octave_NaN;
-
-  for (int i = 0; i < n; i++)
-    {
-      X[i] = args(i).array_value ();
-      Y[i] = args(n+i+1).array_value ();
-
-      if (error_state)
-	{
-	  print_usage ();
-	  return retval;
-	}
-
-      y[i] = Y[i].data ();
-      size[i] =  V.dims()(i);
-
-      if (Y[0].dims () != Y[i].dims ())
-	{
-	  error ("interpn: incompatible size of argument number %d", n+i+2);
-	  return retval;
-	}
-    }
-
-  // offset in memory of each dimension
-
-  scale[0] = 1;
-
-  for (int i = 1; i < n; i++)
-    scale[i] = scale[i-1] * size[i-1];
-
-  // tests if X[0] is a vector, if yes, assume that all elements of X are
-  // in the ndgrid format.
-
-  if (! isvector (X[0]))
-    {
-      for (int i = 0; i < n; i++)
-	{
-	  if (X[i].dims () != V.dims ())
-	    {
-	      error ("interpn: incompatible size of argument number %d", i+1);
-	      return retval;
-	    }
-	  else
-	    {
-              NDArray tmp = NDArray (dim_vector (size[i], 1));
-
-	      for (octave_idx_type j = 0; j < size[i]; j++)
-		tmp(j) =  X[i](scale[i]*j);
-
-              X[i] = tmp;
-	    }
-	}
-    }
-
-  for (int i = 0; i < n; i++)
-    {
-      if (! isvector (X[i]) && X[i].numel () != size[i])
-	{
-	  error ("interpn: incompatible size of argument number %d", i+1);
-	  return retval;
-	}
-      else
-	x[i] = X[i].data ();
-    }
-
-  lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi);
-
-  retval = Vi;
-
-  return retval;
-}

--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -52,12 +52,13 @@
 	dassl.cc det.cc dispatch.cc eig.cc expm.cc fft.cc fft2.cc \
 	fftn.cc fftw.cc filter.cc find.cc fsolve.cc \
 	gammainc.cc gcd.cc getgrent.cc getpwent.cc getrusage.cc \
-	givens.cc hess.cc interpn.cc inv.cc kron.cc lpsolve.cc lsode.cc \
+	givens.cc hess.cc inv.cc kron.cc lpsolve.cc lsode.cc \
 	lu.cc luinc.cc matrix_type.cc md5sum.cc minmax.cc pinv.cc qr.cc \
 	quad.cc qz.cc rand.cc regexp.cc schur.cc sort.cc sparse.cc \
 	spchol.cc spdet.cc spfind.cc spkron.cc splu.cc spparms.cc spqr.cc \
 	sqrtm.cc svd.cc syl.cc symrcm.cc time.cc urlwrite.cc __contourc__.cc \
-	__gnuplot_raw__.l __glpk__.cc __pchip_deriv__.cc __qp__.cc
+	__gnuplot_raw__.l __glpk__.cc __lin_interpn__.cc __pchip_deriv__.cc \
+	__qp__.cc
 
 DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))