function Y = polygen(X, m)
% Y = polygen(X, m)
% generates a matrix of polynomial products of 
% order m from the data matrix X, consisting of:
%    col 1: x1
%    col 2: x2
%    col 3: y
% x1 and x2 are the spatial coordinates, y the 
% observed variable at that coordinate.
% Written by Gerry Middleton, Feb 1997
[r,c] = size(X);
cmax = 0;            % calculate no of terms
for i = 2:m+1 
   cmax = cmax+i;       
end  
Y = zeros(r, cmax);  % reserve memory space
A = zeros(r, cmax);
X1 = zeros(r,cmax);
Y = X(:,1:2);        % starting values
A = X(:,1:2);      % A will be new cols to be added to Y
co = 2;
cn = 2;
for i = 2:m
   X1 = X(:,1)*ones(1,co);  % matrix with cols of x1
   A = X1 .* A;      % elements of A times elements of x1
   A = [A, Y(:,cn).* X(:,2)]; % add column of x2^m
   Y = [Y, A];       % assembled matrix 
   co = co+1;
   cn = cn+co;
end
   
