function b = trend3(X,m,lab,ngrid,reg)
% b = trend3(X,m,lab,ngrid,reg)
% determine the trend surface of order m for a data matrix X 
% with N rows and 3 columns: values of x1, x2, and y
% also plots contours on the original surface, the trend
% surface, and the residuals; and a mesh of the trend. 
% For labels on contours, set lab to 1 (auto) or 2 (manual)
% ngrid is number of grid intervals,e.g, ngrid set to 10
% generates 11 x 11 grid, which is the default
% set reg to 1 if data on a regular grid, otherwise 0 (default)
% NOTE ** coordinates of regular grid must be integers > 0
% If the lowest integer is not 1, the coordinates are rescaled
% Uses polygen, cotomat, dattogr, grtodat
% Written by Gerry Middleton, Feb 1997
if nargin < 5, reg = 0; end
if nargin < 4, ngrid = 10; end
if nargin < 3, lab = 0; end
[N,nc] = size(X);  % determine no of data (rows)
X1 = polygen(X(:,1:2), m);  % generate polynomial terms
X1 = [X1 X(:,3)];  % add dependent variable 
[b,SR,t] = mreg2(X1);
X1 = [ones(N,1) X1];  % expand by adding col of ones
[cols,n] = size(X1);
yhat = X1(:,1:n-1)*b;  % col. vector of estimated values
e = X1(:,n) - yhat;    % col. vector of residuals
if reg == 0
   x1min = ceil(min(X(:,1)));   % find limits of data
   x1max = floor(max(X(:,1)));  % limits are just inside data
   x2min = ceil(min(X(:,2)));   % to help curtail edge effects
   x2max = floor(max(X(:,2)));
   inc = 1/ngrid;             % generate grid coordinates   
   inc1 = (x1max - x1min)*inc;
   inc2 = (x2max - x2min)*inc;
   t1 = x1min:inc1:x1max;
   t2 = x2min:inc2:x2max;
   ng1 = ngrid+1;
   ng2 = ng1;
   [Xi Yi] = meshgrid(t1,t2);
   Z = griddata(X(:,1), X(:,2), X(:,3), Xi, Yi);
else
   eps = 0.000001;
   x1min = min(X(:,1));  % determine max and min for x1
   x1max = max(X(:,1));
   x2min = min(X(:,2));  % and for x2
   x2max = max(X(:,2));
   xx = sort(X(:,1));    % determine increments
   inc1 = max(diff(xx))
   xx = sort(X(:,2));
   inc2 = max(diff(xx))
   t1 = x1min:inc1:x1max+eps; % make x1 and x2 vectors
   t2 = x2min:inc2:x2max+eps; % use eps to avoid rounding errors
   ng1 = length(t1);      % find their size; ng1 is no of cols
   ng2 = length(t2);      % ng2 is no of rows
   ci = round((X(:,1)-x1min)/inc1 + 1); % make integer indices
   ri = round((X(:,2)-x2min)/inc2 + 1);
   Z = full(sparse(ri,ci,X(:,3))); % make array
end
zmin = min(X(:,3));
zmax = max(X(:,3));
v = contvec(zmin,zmax);
vc = length(v);
fprintf('Max value in data: %4.1f  ',zmax);
fprintf('Min value  %4.1f  \n',zmin);
fprintf('There are %2.0f equally spaced contour lines \n',vc);
fprintf('Press <Enter> to continue, at each pause \n');
fprintf('\n');
pause; 
c = contour(Z,v);
if lab == 1 
   clabel(c);
elseif lab == 2
   clabel(c,'manual')
end
xlabel('x'), ylabel('y');
axis image; % change to axis('equal') for MATLAB 4
title('Contours of Original Data');
pause;
if reg == 1
   figure;
   mesh(Z);
   title('Mesh diagram of original data');
   xlabel('X1'),ylabel('X2');
   pause;
end
yhmax = max(yhat);
yhmin = min(yhat);
figure;
fprintf('Max value in trend (at original pts): %4.1f  ',yhmax);
fprintf('Min value  %4.1f  \n',yhmin);
fprintf('\n');
T = cotomat(t1,t2);  % generate list of coords
Xt = polygen(T, m);  % generate polynomial terms
[rt ct] = size(Xt);
Xt = [ones(rt,1) Xt];
that = Xt*b;         % trend at each grid point
Y = dattogr([T that], ng1, ng2);
v = contvec(yhmin,yhmax);  
c = contour(Y,v);
if lab == 1
   clabel(c);
elseif lab == 2
   clabel(c,'manual')
end
axis image; % change to axis('equal') for MATLAB 4
s = num2str(m);
title(['Contours on order ',s,' Trend']);
xlabel('x'), ylabel('y');
pause;
figure;
mesh(Y);
xlabel('X1'),ylabel('X2');
pause;
figure; 
E = Z - Y;
emax = max(max(E));
emin = min(min(E));
v = contvec(emin,emax);
c = contour(E,v);
if lab == 1
   clabel(c);
elseif lab == 2
   clabel(c,'manual');
end
axis image; % change to axis('equal') for MATLAB 4
title('Contours on grid residuals');
xlabel('x'), ylabel('y');
fprintf('Max value in residuals (at original pts): %4.1f  ', max(e));
fprintf('Min value  %4.1f  \n',min(e));
fprintf('\n');
