function b = trend(X,m,lab,ngrid)
% b = trend(X,m,lab,ngrid)
% 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
% Set lab =1 for labels on contours
% ngrid is number of grid intervals, e.g., ngrid = 10
% generates 11x11 node grid, which is the default
% Uses polygen, cotomat, dattogr
% Written by Gerry Middleton, October 1996
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 = [ones(N,1) X1 X(:,3)]; % expand by adding a first column of ones
S = X1'*X1;          % calculate sums of squares and products matrix
[n,n2] = size(S);
SS = S(1:n-1,1:n-1);   % normal equations are  SS*b = y 
y = S(1:n-1,n);
b = SS\y;              % solution of normal equations
yhat = X1(:,1:n-1)*b;  % col. vector of expected values (y hat)
e = X1(:,n) - yhat;    % col. vector of residuals
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;
[Xi Yi] = meshgrid(t1,t2);
Z = griddata(X(:,1), X(:,2), X(:,3), Xi, Yi);
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);
axis image; % change to axis('equal') for MATLAB 5
if lab == 1 
   clabel(c);
end
xlabel('x'), ylabel('y');
title('Contours of Original Data');
pause;
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], ngrid+1, ngrid+1);
v = contvec(yhmin,yhmax);  
c = contour(Y,v);
if lab == 1
   clabel(c);
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; 
E = Z - Y;
emax = max(max(E));
emin = min(min(E));
v = contvec(emin,emax);
c = contour(E,v);
if lab == 1
   clabel(c);
end
axis image; % change to axis('equal') for MATLAB 4
title('Contours on grid residuals');
xlabel('x'), ylabel('y');
pause;
fprintf('Max value in residuals (at original pts): %4.1f  ', max(e));
fprintf('Min value  %4.1f  \n',min(e));
fprintf('\n');
