--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/optiminterp/inst/example_optiminterp.m	Tue Dec 12 23:30:36 2006 +0000
@@ -0,0 +1,39 @@
+% Example program of the optimal interpolation toolbox
+
+
+% number of observations to interpolate
+
+on = 200;
+
+% create randomly located observations within 
+% the square [0 1] x [0 1]
+
+x = rand(1,on);
+y = rand(1,on);
+
+% the underlying function to interpolate
+
+f = sin(6*x) .* cos(6*y);
+
+% the error variance of the observations
+
+var = 0.0 * ones(on,1);
+
+% the grid onto which the observations are interpolated
+
+[xi,yi] = ndgrid(linspace(0,1,100));
+
+% the correlation length in x and y direction
+
+lenx = 0.1;
+leny = 0.1;
+
+% number of influential observations
+
+m = 30;
+
+% run the optimal interpolation
+% fi is the interpolated field and vari is its error variance
+
+[fi,vari] = optiminterp2(x,y,f,var,lenx,leny,m,xi,yi);
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/optiminterp/inst/optiminterp1.m	Tue Dec 12 23:30:36 2006 +0000
@@ -0,0 +1,64 @@
+## Copyright (C) 2006 Alexander Barth
+##
+## This program 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 of the License, or
+## (at your option) any later version.
+##
+## This program 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 this program; if not, write to the Free Software
+## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+## -*- texinfo -*-
+## @deftypefn {Loadable Function} {[@var{fi},@var{vari}] = } optiminterp(@var{x},@var{f},@var{var},@var{lenx},@var{m},@var{xi})
+## Performs a 1D-optimal interpolation (objective analysis).
+##
+## Every elements in @var{f} corresponds to a data point (observation)
+## at location @var{x},@var{y} with the error variance @var{var}.
+##
+## @var{lenx} is correlation length in x-direction.
+## @var{m} represents the number of influential points.
+##
+## @var{xi} is the data points where the field is
+## interpolated. @var{fi} is the interpolated field and @var{vari} is 
+## its error variance.
+##
+## The background field of the optimal interpolation is zero.
+## For a different background field, the background field
+## must be substracted from the observation, the difference 
+## is mapped by OI onto the background grid and finally the
+## background is added back to the interpolated field.
+## @end deftypefn
+
+
+
+function [fi,vari] = optiminterp1(x,f,var,lenx,m,xi)
+
+if (isscalar(var))
+  var = var*ones(size(x));
+end
+
+if isvector(f) & size(f,1) == 1
+   f = f';
+end
+
+if (size(f,1) ~= numel(x)  & numel(f) ~= numel(var))
+  error('optiminterp2: x,var must have the same number of elements');
+end
+
+%whos ox f var lenx m
+
+gx(:,1) = xi(:);
+ox(:,1) = x(:);
+
+%whos ox f var lenx m
+[fi,vari] = optiminterp(ox,f,var,1/lenx,m,gx);
+
+
+%fi = reshape(fi,[numel(xi) size(f,2)]);
+%vari = reshape(vari,size(xi));

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/optiminterp/inst/optiminterp2.m	Tue Dec 12 23:30:36 2006 +0000
@@ -0,0 +1,79 @@
+## Copyright (C) 2006 Alexander Barth
+##
+## This program 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 of the License, or
+## (at your option) any later version.
+##
+## This program 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 this program; if not, write to the Free Software
+## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+## -*- texinfo -*-
+## @deftypefn {Loadable Function} {[@var{fi},@var{vari}] = } optiminterp(@var{x},@var{y},@var{f},@var{var},@var{lenx},@var{leny},@var{m},@var{xi},@var{yi})
+## Performs a 2D-optimal interpolation (objective analysis).
+##
+## Every elements in @var{f} corresponds to a data point (observation)
+## at location @var{x},@var{y} with the error variance @var{var}.
+##
+## @var{lenx} and @var{leny} are correlation length in x-direction
+## and y-direction respectively. 
+## @var{m} represents the number of influential points.
+##
+## @var{xi} and @var{yi} are the data points where the field is
+## interpolated. @var{fi} is the interpolated field and @var{vari} is 
+## its error variance.
+##
+## The background field of the optimal interpolation is zero.
+## For a different background field, the background field
+## must be substracted from the observation, the difference 
+## is mapped by OI onto the background grid and finally the
+## background is added back to the interpolated field.
+## The error variance of the background field is assumed to 
+## have a error variance of one.
+## @end deftypefn
+
+## Copyright (C) 2006, Alexander Barth
+## Author: Alexander Barth <abarth@marine.usf.edu>
+
+function [fi,vari] = optiminterp2(x,y,f,var,lenx,leny,m,xi,yi)
+
+if (isscalar(var))
+  var = var*ones(size(x));
+end
+
+if isvector(f) & size(f,1) == 1
+   f = f';
+end
+
+on = numel(x);
+nf = size(f,3);
+
+if (on*nf ~= numel(f) & on ~= numel(y)  & on ~= numel(var))
+  error('optiminterp2: x,y,f,var must have the same number of elements');
+  return
+end
+
+if (numel(xi) ~= numel(yi))
+  error('optiminterp2: xi and yi must have the same number of elements');
+end
+
+gx(:,1) = xi(:); 
+gx(:,2) = yi(:); 
+
+ox(:,1) = x(:);
+ox(:,2) = y(:);
+
+f=reshape(f,[on nf]);
+
+%whos ox gx f
+
+[fi,vari] = optiminterp(ox,f,var,1./[lenx leny],m,gx);
+
+fi = reshape(fi,[size(xi) nf]);
+vari = reshape(vari,size(xi));

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/optiminterp/inst/optiminterp3.m	Tue Dec 12 23:30:36 2006 +0000
@@ -0,0 +1,82 @@
+## Copyright (C) 2006 Alexander Barth
+##
+## This program 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 of the License, or
+## (at your option) any later version.
+##
+## This program 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 this program; if not, write to the Free Software
+## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+## -*- texinfo -*-
+## @deftypefn {Loadable Function} {[@var{fi},@var{vari}] = } optiminterp(@var{x},@var{y},@var{z},@var{f},@var{var},@var{lenx},@var{leny},@var{lenz},@var{m},@var{xi},@var{yi},@var{zi})
+## Performs a local 3D-optimal interpolation (objective analysis).
+##
+## Every elements in @var{f} corresponds to a data point (observation)
+## at location  @var{x}, @var{y}, @var{z} with the error variance var
+##
+## @var{lenx},@var{leny} and @var{lenz} are correlation length in x-,y- and z-direction
+## respectively. 
+## @var{m} represents the number of influential points.
+##
+## @var{xi},@var{yi} and @var{zi} are the data points where the field is
+## interpolated. @var{fi} is the interpolated field and @var{vari} is 
+## its error variance.
+##
+##
+## The background field of the optimal interpolation is zero.
+## For a different background field, the background field
+## must be substracted from the observation, the difference 
+## is mapped by OI onto the background grid and finally the
+## background is added back to the interpolated field.
+##
+## The error variance of the background field is assumed to 
+## have a error variance of one. 
+## @end deftypefn
+
+## Copyright (C) 2006, Alexander Barth
+## Author: Alexander Barth <abarth@marine.usf.edu>
+
+function [fi,vari] = optiminterp3(x,y,z,f,var,lenx,leny,lenz,m,xi,yi,zi)
+
+if (isscalar(var))
+  var = var*ones(size(x));
+end
+
+if isvector(f) & size(f,1) == 1
+   f = f';
+end
+
+on = numel(x);
+nf = size(f,4);
+
+if (on*nf == numel(f) & on ~= numel(y) & ...
+    on ~= numel(z)  & on ~= numel(var))
+  error('optiminterp3: x,y,z,f,var must have the same number of elements');
+end
+
+if (numel(xi) ~= numel(yi) & numel(xi) ~= numel(zi))
+  error('optiminterp3: xi and yi must have the same number of elements');
+end
+
+gx(:,1) = xi(:); 
+gx(:,2) = yi(:); 
+gx(:,3) = zi(:); 
+
+ox(:,1) = x(:);
+ox(:,2) = y(:);
+ox(:,3) = z(:);
+
+f=reshape(f,[on nf]);
+
+[fi,vari] = optiminterp(ox,f,var,1./[lenx leny lenz],m,gx);
+
+fi = reshape(fi,[size(xi) nf]);
+vari = reshape(vari,size(xi));
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/optiminterp/inst/test_optiminterp.m	Tue Dec 12 23:30:36 2006 +0000
@@ -0,0 +1,112 @@
+% Tests 1D, 2D and 3D optimal interpolation.
+% All tests should pass; any error indicates that either 
+% there is a bug in the optimal interpolation package or 
+% that it is incrorrectly installed.
+
+function test_optiminterp
+
+more off
+
+printf('Testing 1D-optimal interpolation: ');
+
+try
+  % grid of background field
+  xi = linspace(0,1,50);
+  fi_ref = sin( xi*6 );
+
+  % grid of observations
+  x = linspace(0,1,20);
+
+  on = numel(x);
+  var = 0.01 * ones(on,1);
+  f = sin( x*6 );
+
+  m = 15;
+
+  [fi,vari] = optiminterp1(x,f,var,0.1,m,xi);
+
+
+  rms = sqrt(mean((fi_ref(:) - fi(:)).^2));
+
+  if (rms > 0.005) 
+    error('unexpected large difference with reference field');
+  end
+
+  disp('OK');
+
+catch
+  disp('failed');
+  disp(lasterr);
+end
+
+
+printf('Testing 2D-optimal interpolation: ');
+
+try
+  % grid of background field
+  [xi,yi] = ndgrid(linspace(0,1,30));
+  fi_ref = sin( xi*6 ) .* cos( yi*6);
+
+  % grid of observations
+  [x,y] = ndgrid(linspace(0,1,20));
+  x = x(:);
+  y = y(:);
+
+  on = numel(x);
+  var = 0.01 * ones(on,1);
+  f = sin( x*6 ) .* cos( y*6);
+
+  m = 30;
+
+  [fi,vari] = optiminterp2(x,y,f,var,0.1,0.1,m,xi,yi);
+
+
+  rms = sqrt(mean((fi_ref(:) - fi(:)).^2));
+
+  if (rms > 0.005) 
+    error('unexpected large difference with reference field');
+  end
+
+  disp('OK');
+
+catch
+  disp('failed');
+  disp(lasterr);
+end
+
+
+printf('Testing 3D-optimal interpolation: ');
+
+try
+  % grid of background field
+  [xi,yi,zi] = ndgrid(linspace(0,1,15));
+  fi_ref = sin(6*xi) .* cos(6*yi) .* sin(6*zi);
+
+  % grid of observations
+  [x,y,z] = ndgrid(linspace(0,1,10));
+  x = x(:);
+  y = y(:);
+  z = z(:);
+
+  on = numel(x);
+  var = 0.01 * ones(on,1);
+  f = sin(6*x) .* cos(6*y) .* sin(6*z);
+
+  m = 20;
+
+  [fi,vari] = optiminterp3(x,y,z,f,var,0.1,0.1,0.1,m,xi,yi,zi);
+
+
+  rms = sqrt(mean((fi_ref(:) - fi(:)).^2));
+
+  if (rms > 0.04) 
+    error('unexpected large difference with reference field');
+  end
+
+  disp('OK');
+
+catch
+  disp('failed');
+  disp(lasterr);
+end
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/optiminterp/inst/test_optiminterp2.m	Tue Dec 12 23:30:36 2006 +0000
@@ -0,0 +1,119 @@
+% Tests 1D, 2D and 3D optimal interpolation.
+% All tests should pass; any error indicates that either 
+% there is a bug in the optimal interpolation package or 
+% that it is incrorrectly installed.
+
+function test_optiminterp
+
+more off
+
+printf('Testing multiple 1D-optimal interpolation: ');
+
+try
+  % grid of background field
+  xi = linspace(0,1,50)';
+  fi_ref(:,1) = sin( xi*6 );
+  fi_ref(:,2) = cos( xi*6 );
+
+  % grid of observations
+  x = linspace(0,1,20)';
+
+  on = numel(x);
+  var = 0.01 * ones(on,1);
+  f(:,1) = sin( x*6 );
+  f(:,2) = cos( x*6 );
+
+  m = 15;
+
+  [fi,vari] = optiminterp1(x,f,var,0.1,m,xi);
+
+
+  rms = sqrt(mean((fi_ref(:) - fi(:)).^2));
+
+  if (rms > 0.005) 
+    error('unexpected large difference with reference field');
+  end
+
+  disp('OK');
+
+catch
+  disp('failed');
+  disp(lasterr);
+end
+
+
+
+printf('Testing multiple 2D-optimal interpolation: ');
+
+try
+  clear fi_ref f
+  % grid of background field
+  [xi,yi] = ndgrid(linspace(0,1,30));
+
+  fi_ref(:,:,1) = sin( xi*6 ) .* cos( yi*6);
+  fi_ref(:,:,2) = cos( xi*6 ) .* sin( yi*6);
+
+  % grid of observations
+  [x,y] = ndgrid(linspace(0,1,20));
+
+  on = numel(x);
+  var = 0.01 * ones(on,1);
+  f(:,:,1) = sin( x*6 ) .* cos( y*6);
+  f(:,:,2) = cos( x*6 ) .* sin( y*6);
+
+  m = 30;
+
+  [fi,vari] = optiminterp2(x,y,f,var,0.1,0.1,m,xi,yi);
+
+
+  rms = sqrt(mean((fi_ref(:) - fi(:)).^2));
+
+  if (rms > 0.005) 
+    error('unexpected large difference with reference field');
+  end
+
+  disp('OK');
+
+catch
+  disp('failed');
+  disp(lasterr);
+end
+
+printf('Testing multiple 3D-optimal interpolation: ');
+
+
+try
+   clear fi_ref f
+
+  % grid of background field
+  [xi,yi,zi] = ndgrid(linspace(0,1,15));
+  
+  fi_ref(:,:,:,1) = sin(6*xi) .* cos(6*yi) .* sin(6*zi);
+  fi_ref(:,:,:,2) = sin(6*xi) .* cos(6*yi) .* sin(6*zi);
+
+  % grid of observations
+  [x,y,z] = ndgrid(linspace(0,1,10));
+
+  on = numel(x);
+  var = 0.01 * ones(on,1);
+  f(:,:,:,1) = sin(6*x) .* cos(6*y) .* sin(6*z);
+  f(:,:,:,2) = sin(6*x) .* cos(6*y) .* sin(6*z);
+
+  m = 20;
+
+  [fi,vari] = optiminterp3(x,y,z,f,var,0.1,0.1,0.1,m,xi,yi,zi);
+
+
+  rms = sqrt(mean((fi_ref(:) - fi(:)).^2));
+
+  if (rms > 0.04) 
+    error('unexpected large difference with reference field');
+  end
+
+  disp('OK');
+
+catch
+  disp('failed');
+  disp(lasterr);
+end
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/optiminterp/src/Makeconf.in	Tue Dec 12 23:30:36 2006 +0000
@@ -0,0 +1,63 @@
+
+## Makeconf is automatically generated from Makeconf.base and Makeconf.add
+## in the various subdirectories.  To regenerate, use ./autogen.sh to
+## create a new ./Makeconf.in, then use ./configure to generate a new
+## Makeconf.
+
+OCTAVE_FORGE = 1
+
+SHELL = @SHELL@
+
+canonical_host_type = @canonical_host_type@
+prefix = @prefix@
+exec_prefix = @exec_prefix@
+bindir = @bindir@
+mandir = @mandir@
+libdir = @libdir@
+datadir = @datadir@
+infodir = @infodir@
+includedir = @includedir@
+INSTALL = @INSTALL@
+INSTALL_PROGRAM = @INSTALL_PROGRAM@
+INSTALL_SCRIPT = @INSTALL_SCRIPT@
+INSTALL_DATA = @INSTALL_DATA@
+INSTALLOCT=octinst.sh
+
+DESTDIR =
+
+RANLIB = @RANLIB@
+STRIP = @STRIP@
+LN_S = @LN_S@
+
+AWK = @AWK@
+
+# Most octave programs will be compiled with $(MKOCTFILE).  Those which
+# cannot use mkoctfile directly can request the flags that mkoctfile 
+# would use as follows:
+#    FLAG = $(shell $(MKOCTFILE) -p FLAG)
+# The following flags are for compiling programs that are independent
+# of Octave.  How confusing.
+CC = @CC@
+CFLAGS = @CFLAGS@
+CPPFLAGS = @CPPFLAGS@
+CPICFLAG = @CPICFLAG@
+CXX = @CXX@
+CXXFLAGS = @CXXFLAGS@
+CXXPICFLAG = @CXXPICFLAG@
+F77 = @F77@
+FFLAGS = @FFLAGS@
+FPICFLAG = @FPICFLAG@
+
+OCTAVE = @OCTAVE@
+OCTAVE_VERSION = @OCTAVE_VERSION@
+MKOCTFILE = @MKOCTFILE@ -DHAVE_OCTAVE_$(ver) -v
+SHLEXT = @SHLEXT@
+
+MKOCTFILE_FORTRAN_90=@MKOCTFILE_FORTRAN_90@
+OPTIMINTERP_LIBS=@OPTIMINTERP_LIBS@
+OPTIMINTERP_CFLAGS=@OPTIMINTERP_CFLAGS@
+
+%.o: %.c ; $(MKOCTFILE) -c $<
+%.o: %.f ; $(MKOCTFILE) -c $<
+%.o: %.cc ; $(MKOCTFILE) -c $<
+%.oct: %.cc ; $(MKOCTFILE) $<

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/optiminterp/src/configure.base	Tue Dec 12 23:30:36 2006 +0000
@@ -0,0 +1,368 @@
+dnl The configure script is generated by autogen.sh from configure.base 
+dnl and the various configure.add files in the source tree.  Edit 
+dnl configure.base and reprocess rather than modifying ./configure.
+
+dnl autoconf 2.13 certainly doesn't work! What is the minimum requirement?
+AC_PREREQ(2.2)
+
+AC_INIT(configure.base)
+
+PACKAGE=octave-forge
+MAJOR_VERSION=0
+MINOR_VERSION=1
+PATCH_LEVEL=0
+
+dnl Kill caching --- this ought to be the default
+define([AC_CACHE_LOAD], )dnl
+define([AC_CACHE_SAVE], )dnl
+
+dnl uncomment to put support files in another directory
+dnl AC_CONFIG_AUX_DIR(admin)
+
+VERSION=$MAJOR_VERSION.$MINOR_VERSION.$PATCH_LEVEL
+AC_SUBST(PACKAGE)
+AC_SUBST(VERSION)
+
+dnl need to find admin files, so keep track of the top dir.
+TOPDIR=`pwd`
+AC_SUBST(TOPDIR)
+
+dnl if mkoctfile doesn't work, then we need the following:
+dnl AC_PROG_CXX
+dnl AC_PROG_F77
+
+dnl Need C compiler regardless so define it in a way that
+dnl makes autoconf happy and we can override whatever we
+dnl need with mkoctfile -p.
+dnl XXX FIXME XXX should use mkoctfile to get CC and CFLAGS
+AC_PROG_CC
+
+dnl XXX FIXME XXX need tests for -p -c -s in mkoctfile.
+
+dnl *******************************************************************
+dnl Sort out mkoctfile version number and install paths
+
+dnl XXX FIXME XXX latest octave has octave-config so we don't
+dnl need to discover things here.  Doesn't have --exe-site-dir
+dnl but defines --oct-site-dir and --m-site-dir
+
+dnl Check for mkoctfile
+AC_CHECK_PROG(MKOCTFILE,mkoctfile,mkoctfile)
+test -z "$MKOCTFILE" &&	AC_MSG_WARN([no mkoctfile found on path])
+
+AC_SUBST(ver)
+AC_SUBST(subver)
+AC_SUBST(mpath)
+AC_SUBST(opath)
+AC_SUBST(xpath)
+AC_SUBST(altpath)
+AC_SUBST(altmpath)
+AC_SUBST(altopath)
+
+AC_ARG_WITH(path, 
+	[  --with-path             install path prefix],
+	[ path=$withval ])
+AC_ARG_WITH(mpath,
+	[  --with-mpath            override path for m-files],
+	[mpath=$withval])
+AC_ARG_WITH(opath,
+	[  --with-opath            override path for oct-files],
+	[opath=$withval])
+AC_ARG_WITH(xpath,
+	[  --with-xpath            override path for executables],
+	[xpath=$withval])
+AC_ARG_WITH(altpath, 
+	[  --with-altpath          alternative functions install path prefix],
+	[ altpath=$withval ])
+AC_ARG_WITH(altmpath,
+	[  --with-altmpath         override path for alternative m-files],
+	[altmpath=$withval])
+AC_ARG_WITH(altopath,
+	[  --with-altopath         override path for alternative oct-files],
+	[altopath=$withval])	
+
+if test -n "$path" ; then
+   test -z "$mpath" && mpath=$path 
+   test -z "$opath" && opath=$path/oct 
+   test -z "$xpath" && xpath=$path/bin
+   test -z "$altpath" && altpath=$path-alternatives
+fi
+
+if test -n "$altpath" ; then
+   test -z "$altmpath" && altmpath=$altpath 
+   test -z "$altopath" && altopath=$altpath/oct 
+fi
+
+dnl Don't query if path/ver are given in the configure environment
+#if test -z "$mpath" || test -z "$opath" || test -z "$xpath" || test -z "$altmpath" || test -z "$altopath" || test -z "$ver" ; then
+if test -z "$mpath" || test -z "$opath" || test -z "$xpath" || test -z "$ver" ; then
+   dnl Construct program to get mkoctfile version and local install paths
+   cat > conftest.cc <<EOF
+#include <octave/config.h>
+#include <octave/version.h>
+#include <octave/defaults.h>
+
+#define INFOV "\nINFOV=" OCTAVE_VERSION "\n"
+
+#define INFOH "\nINFOH=" OCTAVE_CANONICAL_HOST_TYPE "\n"
+
+#ifdef OCTAVE_LOCALVERFCNFILEDIR
+# define INFOM "\nINFOM=" OCTAVE_LOCALVERFCNFILEDIR "\n"
+#else
+# define INFOM "\nINFOM=" OCTAVE_LOCALFCNFILEPATH "\n"
+#endif
+
+#ifdef OCTAVE_LOCALVEROCTFILEDIR
+# define INFOO "\nINFOO=" OCTAVE_LOCALVEROCTFILEDIR "\n"
+#else
+# define INFOO "\nINFOO=" OCTAVE_LOCALOCTFILEPATH  "\n"
+#endif
+
+#ifdef OCTAVE_LOCALVERARCHLIBDIR
+# define INFOX "\nINFOX=" OCTAVE_LOCALVERARCHLIBDIR  "\n"
+#else
+# define INFOX "\nINFOX=" OCTAVE_LOCALARCHLIBDIR  "\n"
+#endif
+
+const char *infom = INFOM;
+const char *infoo = INFOO;
+const char *infox = INFOX;
+const char *infoh = INFOH;
+const char *infov = INFOV;
+EOF
+
+   dnl Compile program perhaps with a special version of mkoctfile
+   $MKOCTFILE conftest.cc || AC_MSG_ERROR(Could not run $MKOCTFILE)
+
+   dnl Strip the config info from the compiled file
+   eval `strings conftest.o | grep "^INFO.=" | sed -e "s,//.*$,,"`
+   rm -rf conftest*
+
+   dnl set the appropriate variables if they are not already set
+   ver=`echo $INFOV | sed -e "s/\.//" -e "s/\..*$//"`
+   subver=`echo $INFOV | sed -e "[s/^[^.]*[.][^.]*[.]//]"`
+   alt_mbase=`echo $INFOM | sed -e "[s,\/[^\/]*$,,]"`
+   alt_obase=`echo $INFOO | sed -e "[s,/site.*$,/site,]"`
+   test -z "$mpath" && mpath=$INFOM/octave-forge
+   test -z "$opath" && opath=$INFOO/octave-forge
+   test -z "$xpath" && xpath=$INFOX
+   test -z "$altmpath" && altmpath=$alt_mbase/octave-forge-alternatives/m
+   test -z "$altopath" && altopath=$alt_obase/octave-forge-alternatives/oct/$INFOH
+fi
+
+dnl *******************************************************************
+
+dnl XXX FIXME XXX Should we allow the user to override these?
+dnl Do we even need them?  The individual makefiles can call mkoctfile -p
+dnl themselves, so the only reason to keep them is for configure, and
+dnl for those things which are not built using mkoctfile (e.g., aurecord)
+dnl but it is not clear we should be using octave compile flags for those.
+
+dnl C compiler and flags
+AC_MSG_RESULT([retrieving compile and link flags from $MKOCTFILE])
+CC=`$MKOCTFILE -p CC`
+CFLAGS=`$MKOCTFILE -p CFLAGS`
+CPPFLAGS=`$MKOCTFILE -p CPPFLAGS`
+CPICFLAG=`$MKOCTFILE -p CPICFLAG`
+LDFLAGS=`$MKOCTFILE -p LDFLAGS`
+LIBS=`$MKOCTFILE -p LIBS`
+AC_SUBST(CC)
+AC_SUBST(CFLAGS)
+AC_SUBST(CPPFLAGS)
+AC_SUBST(CPICFLAG)
+
+dnl Fortran compiler and flags
+F77=`$MKOCTFILE -p F77`
+FFLAGS=`$MKOCTFILE -p FFLAGS`
+FPICFLAG=`$MKOCTFILE -p FPICFLAG`
+AC_SUBST(F77)
+AC_SUBST(FFLAGS)
+AC_SUBST(FPICFLAG)
+
+dnl C++ compiler and flags
+CXX=`$MKOCTFILE -p CXX`
+CXXFLAGS=`$MKOCTFILE -p CXXFLAGS`
+CXXPICFLAG=`$MKOCTFILE -p CXXPICFLAG`
+AC_SUBST(CXX)
+AC_SUBST(CXXFLAGS)
+AC_SUBST(CXXPICFLAG)
+
+dnl *******************************************************************
+
+dnl Check for features of your version of mkoctfile.
+dnl All checks should be designed so that the default
+dnl action if the tests are not performed is to do whatever
+dnl is appropriate for the most recent version of Octave.
+
+dnl Define the following macro:
+dnl    OF_CHECK_LIB(lib,fn,true,false,helpers)
+dnl This is just like AC_CHECK_LIB, but it doesn't update LIBS
+AC_DEFUN(OF_CHECK_LIB,
+[save_LIBS="$LIBS"
+AC_CHECK_LIB($1,$2,$3,$4,$5)
+LIBS="$save_LIBS"
+])
+
+dnl Define the following macro:
+dnl    TRY_MKOCTFILE(msg,program,action_if_true,action_if_false)
+dnl
+AC_DEFUN(TRY_MKOCTFILE,
+[AC_MSG_CHECKING($1)
+cat > conftest.cc << EOF
+#include <octave/config.h>
+$2
+EOF
+ac_try="$MKOCTFILE -c conftest.cc"
+if AC_TRY_EVAL(ac_try) ; then
+   AC_MSG_RESULT(yes)
+   $3
+else
+   AC_MSG_RESULT(no)
+   $4
+fi
+])
+
+dnl
+dnl Check if F77_FUNC works with MKOCTFILE
+dnl
+TRY_MKOCTFILE([for F77_FUNC],
+[int F77_FUNC (hello, HELLO) (const int &n);],,
+[MKOCTFILE="$MKOCTFILE -DF77_FUNC=F77_FCN"])
+
+dnl **********************************************************
+
+dnl Evaluate an expression in octave
+dnl
+dnl OCTAVE_EVAL(expr,var) -> var=expr
+dnl
+AC_DEFUN(OCTAVE_EVAL,
+[AC_MSG_CHECKING([for $1 in Octave])
+$2=`echo "disp($1)" | $OCTAVE -qf`
+AC_MSG_RESULT($$2)
+AC_SUBST($2)
+])
+
+dnl Check status of an octave variable
+dnl
+dnl OCTAVE_CHECK_EXIST(variable,action_if_true,action_if_false)
+dnl
+AC_DEFUN(OCTAVE_CHECK_EXIST,
+[AC_MSG_CHECKING([for $1 in Octave])
+if test `echo 'disp(exist("$1"))' | $OCTAVE -qf`X != 0X ; then
+   AC_MSG_RESULT(yes)
+   $2
+else
+   AC_MSG_RESULT(no)
+   $3
+fi
+])
+
+dnl should check that $(OCTAVE) --version matches $(MKOCTFILE) --version
+AC_CHECK_PROG(OCTAVE,octave,octave)
+OCTAVE_EVAL(OCTAVE_VERSION,OCTAVE_VERSION)
+
+dnl grab canonical host type so we can write system specific install stuff
+OCTAVE_EVAL(octave_config_info('canonical_host_type'),canonical_host_type)
+
+dnl grab SHLEXT from octave config
+OCTAVE_EVAL(octave_config_info('SHLEXT'),SHLEXT)
+
+AC_PROG_LN_S
+
+AC_PROG_RANLIB
+
+dnl Use $(COPY_FLAGS) to set options for cp when installing .oct files.
+COPY_FLAGS="-Rfp"
+case "$canonical_host_type" in
+    *-*-linux*)
+        COPY_FLAGS="-fdp"
+    ;;
+esac
+AC_SUBST(COPY_FLAGS)
+
+dnl Use $(STRIP) in the makefile to strip executables.  If not found, 
+dnl STRIP expands to ':', which in the makefile does nothing.
+dnl Don't need this for .oct files since mkoctfile handles them directly
+STRIP=${STRIP-strip}
+AC_CHECK_PROG(STRIP,$STRIP,$STRIP,:)
+
+dnl Strip on windows, don't strip on Mac OS/X or IRIX
+dnl For the rest, you can force strip using MKOCTFILE="mkoctfile -s"
+dnl or avoid strip using STRIP=: before ./configure
+case "$canonical_host_type" in
+    powerpc-apple-darwin*|*-sgi-*)
+	STRIP=:
+    ;;
+    *-cygwin-*|*-mingw-*) 
+	MKOCTFILE="$MKOCTFILE -s" 
+    ;;
+esac
+
+AC_SUBST(MKOCTFILE_FORTRAN_90)
+AC_SUBST(OPTIMINTERP_LIBS)
+AC_SUBST(OPTIMINTERP_CFLAGS)
+
+OPTIMINTERP_LIBS="$LDFLAGS"
+OPTIMINTERP_CFLAGS="$CPPFLAGS"
+
+dnl By default, autoconf tries to compile Fortran programs with
+dnl optimization (-O2) which mkoctfile does not support as argument
+
+FCFLAGS=
+
+dnl Set Language to Fortran
+AC_LANG(Fortran)
+
+dnl Use mkoctfile as Fortran 90 compiler.
+dnl No tests are performed at this stage if mkoctfile can actually handle
+dnl Fortran 90 code.
+
+AC_PROG_FC([mkoctfile], 1990)
+
+dnl Test if mkoctfile and the Fortran 90 compiler support
+dnl the extention .F90 (free-form with pre-processor directives)
+
+AC_FC_SRCEXT(F90,,AC_ERROR([mkoctfile does not accept files with the extention .F90. \
+Support for this file extension has been added to version 2.9.9 of octave.]))
+
+dnl Test if mkoctfile can handle Fortran 90 features such as modules
+
+AC_MSG_CHECKING([whether mkoctfile accepts Fortran 90])
+
+AC_COMPILE_IFELSE([[\
+module test_module; \
+contains; \
+end module \
+]], MKOCTFILE_FORTRAN_90=yes, MKOCTFILE_FORTRAN_90=no)
+
+OPTIMINTERPSTATUS=$MKOCTFILE_FORTRAN_90
+AC_MSG_RESULT([$MKOCTFILE_FORTRAN_90])
+
+if test $MKOCTFILE_FORTRAN_90 = no ; then
+    AC_MSG_ERROR([mkoctfile does not accept Fortran 90 code. Octave was not build with a Fortran 90 compiler.])
+fi
+
+dnl Delete the module file. Capitalization is compiler dependent
+rm -f test_module.mod TEST_MODULE.mod test_module.MOD TEST_MODULE.MOD
+
+
+CONFIGURE_OUTPUTS="Makeconf"
+STATUS_MSG="
+octave commands will install into the following directories:
+   m-files:   $mpath
+   oct-files: $opath
+   binaries:  $xpath
+alternatives:
+   m-files:   $altmpath
+   oct-files: $altopath
+
+shell commands will install into the following directories:
+   binaries:  $bindir
+   man pages: $mandir
+   libraries: $libdir
+   headers:   $includedir
+
+octave-forge is configured with
+   octave:      $OCTAVE (version $OCTAVE_VERSION)
+   mkoctfile:	$MKOCTFILE for Octave $subver
+   optiminterp toolbox: $OPTIMINTERPSTATUS"

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/optiminterp/src/example_optiminterp.F90	Tue Dec 12 23:30:36 2006 +0000
@@ -0,0 +1,100 @@
+!  
+! Example for using the n-dimentional optimal interpolation Fortran module
+! Copyright (C) 2005 Alexander Barth <abarth@marine.usf.edu>
+!  
+! This program 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
+! of the License, or (at your option) any later version.
+!  
+! This program 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 this program; if not, write to the Free Software
+! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, 
+! MA 02110-1301, USA.
+!  
+
+! Author: Alexander Barth <abarth@marine.usf.edu>
+
+program test_optiminterp
+ use optimal_interpolation
+ implicit none
+
+ ! n times n is the sized of the 2D background grid
+ ! on: number of observations
+ ! m: number of influential observations
+
+ integer, parameter :: n=100, on=200, m=30
+
+ ! xi(1,:) and xi(2,:) represent the x and y coordindate of the 
+ ! grid of the interpolated field
+ ! fi and vari are the interpolated field and its error variance resp.
+
+ real(wp) :: xi(2,n*n), fi(1,n*n), vari(n*n)
+
+ ! xi(1,:) and xi(2,:) represent the x and y coordindate of the 
+ ! observations
+ ! f and var are observations and their error variance resp.
+
+
+ real(wp) :: x(2,on), var(on), f(1,on)
+
+ ! param: inverse of the correlation length
+
+ real(wp) :: param(2)
+
+ integer :: i,j,l
+
+ ! we use a simple random generator instead of Fortran 90's
+ ! random_number in the hope that the results will be
+ ! platform independent
+
+ integer(8), parameter :: A = 1664525_8, B = 1013904223_8,
+ &              Mo = 4294967296_8
+ integer(8) :: next = 0_8
+
+
+ ! create a regular 2D grid
+ l = 1
+ do j=1,n
+   do i=1,n
+     xi(1,l) = (i-1.)/(n-1.)
+     xi(2,l) = (j-1.)/(n-1.)
+     l = l+1
+   end do
+ end do
+
+ ! param is the inverse of the correlation length
+ param = 1./0.1
+
+ ! the error variance of the observations
+ var = 0.01
+
+ ! location of observations
+
+ do j=1,on
+   do i=1,2
+     ! simple random generator
+     next = mod(A*next  + B,Mo)
+     x(i,j) = real(next,8)/real(Mo,8)
+   end do
+ end do
+
+ ! the underlying function to interpolate
+
+ f(1,:) = sin( x(1,:)*6 ) * cos( x(2,:)*6);
+
+ ! fi is the interpolated function and vari its error variance
+ call optiminterp(x,f,var,param,m,xi,fi,vari)
+
+ ! control value
+
+ write(6,'(A,F10.6)') 'Expected:',2.2205936104348591E-002
+ write(6,'(A,F10.6)') 'Obtained:',fi(1,1)
+
+end program test_optiminterp
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/optiminterp/src/optimal_interpolation.F90	Tue Dec 12 23:30:36 2006 +0000
@@ -0,0 +1,345 @@
+
+!  Fortran 90 module for n-dimensional optimal interpolation.
+!  Copyright (C) 2005 Alexander Barth <abarth@marine.usf.edu>
+!  
+!  This program 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
+!  of the License, or (at your option) any later version.
+!  
+!  This program 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 this program; if not, write to the Free Software
+!  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+!  
+
+!  Author: Alexander Barth <abarth@marine.usf.edu>
+!  Dependencies: LAPACK (dsyev)
+
+!  David Saunders (NASA Ames Research Center)
+!     optimizations and code clean-up
+
+#define DIAG_OBS_COVAR
+      module optimal_interpolation
+!     working precision 
+!     4 = simple precision, 8 is double precision
+
+      integer, parameter :: wp = 8
+
+
+      contains
+
+!     --------------------------------------------------------------------------
+
+      subroutine select_nearest(x,ox,param,m,index,distance)
+
+!     Select the m observations from ox(1:nd,1:on) closest to point x(1:nd).
+
+!     Arguments:
+
+      implicit none
+      real(wp),intent(in)  :: x(:),ox(:,:),param(:)
+      integer, intent(in)  :: m
+      real(wp),intent(out) :: distance(m)
+      integer, intent(out) :: index(m)
+
+!     Local variables:
+
+      real(wp) :: d(size(ox,2))
+      integer  :: i
+
+!     Execution:
+
+!     Calculate a measure of (squared) distance to each observation:
+
+      do i=1,size(ox,2)
+        d(i) = sum(((x - ox(:,i)) * param)**2)
+      end do
+
+      call sort(d,m,index)
+
+      distance = d(index)
+
+      end subroutine select_nearest
+
+
+!     --------------------------------------------------------------------------
+
+      subroutine sort(d,m,pannier)
+
+!     Return the indices of the m smallest elements in d(:).
+!     The algorithm is succinctly coded, but would a heap sort be faster?
+
+!     Arguments:
+
+      implicit none
+      real(wp), intent(in)  :: d(:)
+      integer,  intent(in)  :: m
+      integer,  intent(out) :: pannier(m)
+
+      integer :: i,max_pannier(m)
+
+      do i=1,m
+        pannier(i) = i
+      end do
+
+      max_pannier = maxloc(d(pannier))
+
+      do i=m+1,size(d)
+        if (d(i) .lt. d(pannier(max_pannier(1)))) then
+          pannier(max_pannier(1)) = i
+          max_pannier = maxloc(d(pannier))
+        end if
+      end do
+
+      end subroutine sort
+
+!     --------------------------------------------------------------------------
+
+      subroutine observation_covariance(ovar,index,R)
+      implicit none
+      real(wp),intent(in) ::ovar(:)
+      integer,intent(in) :: index(:)
+#     ifdef DIAG_OBS_COVAR
+      real(wp),    intent (out) :: R(size(index))
+#     else
+      real(wp),    intent (out) :: R(size(index),size(index))
+      integer :: i
+#     endif
+
+
+#     ifdef DIAG_OBS_COVAR
+      R = ovar(index)
+#     else
+      R = 0
+      do i=1,size(index)
+        R(i,i) = ovar(index(i))
+      end do
+#     endif
+
+      end subroutine observation_covariance
+
+!     --------------------------------------------------------------------------
+
+      function  background_covariance(x1,x2,param) result(c)
+      implicit none
+      real(wp),intent(in) :: x1(:),x2(:),param(:)
+      real(wp) :: c
+
+      real(wp) :: d(size(x1))
+
+      d = (x1 - x2)*param
+
+      c = exp(-sum(d**2))
+
+      end function background_covariance
+
+!     --------------------------------------------------------------------------
+
+      subroutine pinv (A, tolerance, work, D)
+
+!     Compute pseudo-inverse A+ of symmetric A in factored form U D+ U', where
+!     U overwrites A and D is diagonal matrix D+.
+
+!     Saunders notes: Working with the factors of the pseudo-inverse is
+!                     preferable to multiplying them together (more stable,
+!                     less arithmetic).  Also, the choice of tolerance is
+!                     not straightforward.  If A is noisy, try 1.e-2 * || A ||,
+!                     else machine-eps. * || A ||  (1 norms).
+!     Arguments:
+
+      real(wp), intent (inout) :: A(:,:)  ! Upper triangle input; orthogonal U out
+      real(wp), intent (in)    :: tolerance
+      real(wp), intent (out)   :: work(:), D(:)
+
+!     Local variables:
+
+      integer :: i, j, k, info, N
+
+!     Execution:
+
+      N = size (A,1)
+
+!     Eigendecomposition/SVD of symmetric A:
+
+      call dsyev ('V', 'U', N, A, N, D, work, size (work), info)
+
+!     Diagonal factor D+ of pseudo-inverse:
+
+      do i = 1, N
+         if (D(i) > tolerance) then
+            D(i) = 1. / D(i)
+         else
+            D(i) = 0.
+         end if
+      end do
+
+
+      end subroutine pinv
+
+!     --------------------------------------------------------------------------
+
+      function pinv_workspace(N) result(lwork)
+
+!     Determine the workspace needed for dsyev.
+
+      implicit none
+      integer,intent(in) :: N
+      integer :: lwork
+
+      integer :: info
+      real(wp) :: dummy,rwork
+
+      call dsyev('V','U', N, dummy,N, dummy, rwork, -1, info )
+      lwork = ceiling(rwork)
+
+      end function
+
+!     --------------------------------------------------------------------------
+
+      subroutine optiminterp(ox,of,ovar,param,m,gx,gf,gvar)
+
+!     Main optimal interpolation routine
+
+      implicit none
+      real(wp), intent(in)  :: ox(:,:), of(:,:)  ! Observations
+      real(wp), intent(in)  :: ovar(:)           ! Observation error variances
+      real(wp), intent(in)  :: param(:)          ! inverse of correlation lengths
+      integer,  intent(in)  :: m                 ! # nearest observations used
+      real(wp), intent(in)  :: gx(:,:)           ! Target grid coords.
+      real(wp), intent(out) :: gf(:,:), gvar(:)  ! Interpolated fields.
+                                                 ! and error variances
+!     Local variables:
+
+      real(wp) :: HPH(m,m), PH(m), iA(m,m), PHiA(m),A(m,m),D(m)
+
+#ifdef DIAG_OBS_COVAR
+      real(wp) :: R(m)
+#else
+      real(wp) :: R(m,m)
+#endif
+
+      integer  :: gn,nf,index(m)
+      real(wp) :: distance(m)
+
+      integer  :: i,j1,j2,j,lwork
+      real(wp) :: tolerance = 1e-5
+
+#     ifdef VERBOSE
+      integer  :: percentage_done
+#     endif
+#     ifdef STATIC_WORKSPACE
+      real(wp) :: work((m+2)*m)
+#     else
+      real(wp), allocatable :: work(:)
+#     endif
+
+!     Execution:
+
+      gn = size(gx,2)
+      nf = size (of, 1)  ! # functions at each observation point
+
+#     ifndef STATIC_WORKSPACE
+!     query and allocate workspace for pseudo-inverse
+      lwork = pinv_workspace(m)
+#     endif
+
+!$omp parallel private(work,i,iA,PHiA,index,distance,HPH,j1,j2)
+
+#     ifndef STATIC_WORKSPACE
+      allocate(work(lwork))
+#     endif
+      
+#     ifdef VERBOSE
+      percentage_done = 0
+#     endif
+
+!$omp do 
+      do i=1,gn
+#       ifdef VERBOSE
+        if (percentage_done .ne. int(100.*real(i)/real(gn))) then
+          percentage_done = int(100.*real(i)/real(gn))
+          write(6,*) 'done ',percentage_done
+        end if
+#       endif
+
+!       get the indexes of the nearest observations
+
+        call select_nearest(gx(:,i),ox,param,m,index,distance)
+     
+!       form compute the error covariance matrix of the observation 
+
+        call observation_covariance(ovar,index,R)
+
+!       form the error covariance matrix background field
+
+        do j2 = 1, m
+          ! Upper triangle only
+
+           do j1 = 1, j2  
+             A(j1,j2) = &
+                 background_covariance(ox(:,index(j1)),ox(:,index(j2)),param)
+           end do
+
+           PH(j2) = background_covariance(gx(:,i),ox(:,index(j2)),param)
+        end do
+
+!       covariance matrix of the innovation
+
+#ifdef DIAG_OBS_COVAR
+        do j = 1, m
+           A(j,j) = A(j,j) + R(j)
+        end do
+#else
+        A  = A + R
+#endif
+
+!       pseudo inverse of the covariance matrix of the innovation
+
+        call pinv(A,tolerance,work,D)
+
+        PHiA = matmul (A, D * matmul (PH, A))
+
+!       compute the analysis for all fields
+
+        gf(:,i) = matmul(of(:,index),PHiA)
+
+!       compute the error variance of the analysis
+
+        gvar(i) = 1. - dot_product(PHiA,PH)
+
+      end do
+!$omp end do
+
+#     ifndef STATIC_WORKSPACE
+      deallocate(work)
+#     endif
+ 
+!$omp end parallel 
+
+
+      end subroutine optiminterp
+
+#ifdef FORTRAN_2003_INTEROP
+
+      subroutine optiminterp_gw(n,gn,on,nparam,ox,of,ovar,
+     &                           param,m,gx,gf,gvar) bind(c)
+      USE ISO_C_BINDING
+      implicit none
+      integer(C_INT) :: m,n,gn,on,nparam
+      real(C_DOUBLE) :: gx(n,gn),ox(n,on),of(on),ovar(on),param(nparam)
+      real(C_DOUBLE) :: gf(gn),gvar(gn)
+
+
+      call optiminterp(ox,of,ovar,param,m,gx,gf,gvar)
+      
+      end subroutine optiminterp_gw
+
+
+#endif
+
+      end module optimal_interpolation

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/optiminterp/src/optiminterp.cc	Tue Dec 12 23:30:36 2006 +0000
@@ -0,0 +1,113 @@
+/*
+  n-dimentional optimal interpolation
+  Copyright (C) 2005 Alexander Barth <abarth@marine.usf.edu>
+
+  This program 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
+  of the License, or (at your option) any later version.
+
+  This program 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 this program; if not, write to the Free Software
+  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+
+
+  Author: Alexander Barth <abarth@marine.usf.edu>
+*/									
+
+
+#include <octave/oct.h>
+#include "f77-fcn.h"
+
+extern "C"
+{
+  int F77_FUNC (optiminterp_wrapper, OPTIMINTERP_WRAPPER) 
+    (
+     const int& n, const int& nf, const int& gn, const int& on, const int& nparam,  
+     double* ox, double* of, double* ovar, double* param,
+     const int& m, double* gx,  double* gf, double* gvar);
+}
+
+DEFUN_DLD (optiminterp, args, ,
+	   "- Loadable Function: [gf,gvar] = optiminterp(ox,of,ovar,param,m,gx)\n\
+     \n\
+     Performs an N-dimensional optimal interpolation.")
+{
+  octave_value_list retval;
+  
+  if (args.length() != 6) {
+    print_usage ();
+    return octave_value();
+  }
+
+  Matrix ox = args(0).matrix_value();
+  Matrix of = args(1).matrix_value();
+  Matrix ovar = args(2).matrix_value();
+  Matrix param = args(3).matrix_value();
+  int m = (int)args(4).scalar_value();
+  Matrix gx = args(5).matrix_value();
+
+  /*
+   In Octave and Matlab convention (see gridatan), the positions of series 
+   of points in n-D space is represented by an m-by-n matrix. However, it is more
+   efficient to represent it as a n-ny-m matrix. The Fortran code uses
+   the later convension. Therefore we need to transpose the following 
+   matrices:
+  */
+
+  of = of.transpose();
+  gx = gx.transpose();
+  ox = ox.transpose();
+  //octave_stdout << "shape(of) " << of.rows() << "x" << of.cols() << std::endl;
+
+  int n = gx.rows();
+  int nf = of.rows();
+  int gn = gx.cols();
+  int on = ox.cols();
+  int nparam = param.numel();
+  Matrix gf = Matrix(nf,gn);
+  Matrix gvar = Matrix(gn,1);
+
+  /*
+  octave_stdout << "nf " << nf << std::endl;
+  octave_stdout << "n "  << n << std::endl;
+  octave_stdout << "gn " << gn << std::endl;
+  octave_stdout << "on " << on << std::endl;
+  */
+
+  OCTAVE_QUIT;
+  if (gx.rows() != ox.rows()) {
+    error ("optiminterp: number of rows in gx and ox must be the same");
+    return octave_value();
+  }
+
+  if (ox.cols() != of.cols() || ox.cols() != ovar.numel()) {
+    error ("optiminterp: number of colums in ox must be equal to the number of elements in of and ovar");
+    return octave_value();
+  }
+
+  if (m > on) {
+    warning("optiminterp: Number of influential points is larger than number of data points. Limiting the number of influential points to the number of data points");
+    m = on;
+  }
+
+  F77_XFCN (optiminterp_wrapper, OPTIMINTERP_WRAPPER,
+	    ( n,nf,gn,on,nparam,ox.fortran_vec(),of.fortran_vec(),ovar.fortran_vec(),
+              param.fortran_vec(),m,gx.fortran_vec(),gf.fortran_vec(),gvar.fortran_vec()));
+  if (f77_exception_encountered) {
+      error ("unrecoverable error in optiminterp");
+      return retval;
+  }
+
+  gf = gf.transpose();
+
+  retval(0) = gf;
+  retval(1) = gvar;
+
+  return retval;
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/optiminterp/src/optiminterp_wrapper.F90	Tue Dec 12 23:30:36 2006 +0000
@@ -0,0 +1,42 @@
+!
+!  n-dimentional optimal interpolation
+!  Copyright (C) 2005 Alexander Barth <abarth@marine.usf.edu>
+!
+!  This program 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
+!  of the License, or (at your option) any later version.
+!
+!  This program 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 this program; if not, write to the Free Software
+!  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+!
+
+!  Author: Alexander Barth <abarth@marine.usf.edu>
+
+subroutine optiminterp_wrapper(n,nf,gn,on,nparam,ox,of,ovar,    &
+                               param,m,gx,gf,gvar)
+ use optimal_interpolation
+ implicit none
+ integer :: m,n,nf,gn,on,nparam
+ real(8) :: gx(n,gn),ox(n,on),of(nf,on),ovar(on),param(nparam)
+ real(8) :: gf(nf,gn),gvar(gn)
+
+
+ call optiminterp(ox,of,ovar,param,m,gx,gf,gvar)
+!!$
+!!$ integer :: i,j
+!!$ real :: x=1
+!!$ do i=1,1000
+!!$   write(6,*) 'i',i
+!!$   do j= 1,1000000
+!!$    x = x + sin(x)
+!!$   end do
+!!$ end do
+
+end subroutine optiminterp_wrapper