## Copyright (C) 2007 Lei Jia
##
## 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} {} interpolate (@var{z}, @var{x}, @var{val})
## This function uses linear interpolation to return a cross-sectional
## set of data from the mesh data.
##
## @var{z} is the mesh data matrix.
##
## @var{x} is the column vector or row vector for the matrix.
##
## @itemize @bullet
##
## @item
## When @var{x} is a column vector, the return value will be
## a row-wise interpolation of the mesh matrix.
##
## @item
## When @var{x} is a row vector, the return value will be
## a column-wise interpolation of the mesh matrix.
## @end itemize
##
## @var{val} is the value to be interpolated.
## @end deftypefn

## Author: Lei Jia

function inter = interpolate(z, x, val)

	row = rows(x);
	col = columns(x);
	
	if (row == 1)
		i = (col == columns(z));
	elseif (col == 1)
		i = (row == rows(z));
	else
		error("x must be either a row vector or a column vector");
	endif
	
	if (row == 1 && i == 1)
		for ic = 1:col
			if ((x(ic) == val) || (ic == 1 && x(ic) > val) || (ic == col && x(ic) < val))
				for ir = 1:rows(z)
					inter(ir) = z(ir, ic);
				endfor
				break;
			elseif (ic != 1 && x(ic) > val && x(ic - 1) < val)
				## interpolation
				for ir = 1:rows(z)
					inter(ir) = z(ir, ic - 1) + (val - x(ic - 1)) * (z(ir, ic) - z(ir, ic - 1)) / (x(ic) - x(ic - 1));
				endfor
				break;
			endif
		endfor
		inter = inter';
	elseif (col == 1 && i == 1)
		for ir = 1:row
			if ((x(ir) == val) || (ir == 1 && x(ir) > val) || (ir == row && x(ir) < val))
				for ic = 1:columns(z)
					inter(ic) = z(ir, ic);
				endfor
				break;
			elseif (ir != 1 && x(ir) > val && x(ir - 1) < val)
				## interpolation
				for ic = 1:columns(z)
					inter(ic) = z(ir - 1, ic) + (val - x(ir - 1)) * (z(ir, ic) - z(ir - 1, ic)) / (x(ir) - x(ir - 1));
				endfor
				break;
			endif
		endfor
	else
		error("the length of x must be equal to either the rows or the columns of z");
	endif

endfunction