--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/interpn.cc	Thu Feb 15 07:52:28 2007 +0000
@@ -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 "dMatrix.cc"
+
+#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;
+}