function Z = polyval2(C,x,y, method)
% polyval2.m  evaluates a 2D polynomial fitting from intfit2.m
%
% USAGE:  Z = polyval2(C,X,Y,'method') 
% ======  returns matrix Z containing evaluations of a least-squares 
%         polynomial in x and y, 
% where coefficients "C" has been based on 'method' = 
%                    'linear' - bilinear linear squares fitting  OR
%                     'cubic'  - bicubic linear squares fitting
% Again non-equally spaced X and Y are permitted but note that
%       you may not plot the polynomial at nonuniformly meshes.
%
% For example, generate a coarse 2D sine curve and a least squares fitting
% ver finer abscissa:
%                  x = 0:10; y = 1:9;  [x y] = meshgrid(x,y) ; 
%                  z = sin(x.*y);
%                 xi = 0:.25:10; yi=2:.5:8 ; [xi yi]=meshgrid(xi,yi);
%                 C  = polyfit2(x,y,z, 'cubic');
%                 zi = polyval2(C, xi,yi, 'cubic');

if nargin<=2
	disp('Insufficient number of options are supplied');
	disp('                       Error in polyval2')
	help polyval2
	return
end

% Check the arguments.
if nargin<4,
  method = 'linear';
end

if ~isstr(method),   %%% Any string
  error('METHOD must be one of the strings: ''linear'',''cubic''.');
end

if size(x)~=size(y),
  error('X must have the same size as Y.');
end

[m n] = size(x); mn = m*n ;
  x = reshape(x,mn,1) ;
  y = reshape(y,mn,1) ;

if method(1)=='l', % Linear interpolation
  Z = C(1) + C(2)*x + C(3)*y + C(4)*x.*y ;

elseif method(1)=='c', % Cubic interpolation
  Z = C(1) + C(2)*x + C(3)*y + C(4)*x.*y + C(5)*x.^2 + C(6)*y.^2 ;
  Z = Z    + C(7)*(x.^2).*y + C(8)*x.*y.^2 + C(9)*x.^3 + C(10)*y.^3 ;
else
  error([method,' is an invalid method.']);
end

% Recover a matrix form
  x = reshape(x,m,n) ;
  y = reshape(y,m,n) ;
  Z = reshape(Z,m,n) ;
