function C = polyfit2(x,y,z, method)  

% polyfit2 is for 2-D data fitting using least squares

%

% USAGE:   C = polyfit2(X,Y,Z, 'method') 

% ======where an output vector C contains the bi-linear or bi-cubic

%       coefficients of a least-squares polynomial in x and y, and

%       input matrices X, Y, Z are for a 2D function z=f(x,y).

%

% Here 'method' can be

%         'linear' - bilinear least squares fitting

%         'cubic'  -  bicubic least squares fitting

%       Non-equally spaced (or even non-monotonic) X and Y are permitted.

%

%       For example, generate a coarse 2D curve and a least squares fitting

%       over finer mesh (meshdom with Matlab 3.5 BUT meshgrid with Matlab 4)

%                  x = 0:10; y = 1:9;  [x y] = meshdom(x,y) ; 

%                            z = sin(x.*y);

%                  xi = 0:.25:10; yi=2:.5:8 ; [xi yi]=meshdom(xi,yi);

%                  C  = polyfit2(x,y,z, 'cubic');

%                  zi = polyval2(C, xi,yi, 'cubic');

if nargin<=2, format compact

   disp('Insufficient number of options are supplied');

   disp('                       Error in polyfit2');help polyfit2

   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(z),

  error('X must have the same dimension as Z.');                      

end

if size(y)~=size(z),

  error('Y must have the same dimension as Z.');                      

end



[m,n] = size(z); mn = m*n  ;  



  x = reshape(x,mn,1) ;

  y = reshape(y,mn,1) ;

  z = reshape(z,mn,1) ;



if method(1)=='l', % Linear least squares 

  A = [ ones(size(x))  x  y  x.*y ] ;

  C  = A \ z ;



elseif method(1)=='c', % Cubic least squares 

  A = [ ones(size(x))  x  y  x.*y  x.^2  y.^2  (x.^2).*y  x.*(y.^2)  x.^3  y.^3];

  C  = A \ z ;



else

  error([method,' is an invalid method.']);

end

 

%% Finally  Output as a row vector

C  = C' ;

  x = reshape(x,m,n) ;

  y = reshape(y,m,n) ;

  z = reshape(z,m,n) ;



