function [tz] = gridcon(X,tx,ty,opt);
% function [tz] = gridcon(X,tx,ty,opt);
% grids and contours the data in X - an Nx3
% matrix with each row consisting of (x,y,z)
% where x and y are coordinates and z is a
% measured property. tx and ty are row vectors
% defining the grid in the x and y direction
% respectively. opt defines interpolation
% options: 1 = linear; 2 = cubic; 3 = nearest;
% 4 = invdist
if nargin < 4, opt = 1; end
xmin = min(tx);
xmax = max(tx);
ymin = min(ty);
ymax = max(ty);
plot(X(:,1), X(:,2),'o');
axis([xmin xmax ymin ymax]);
axis('equal');
hold on;
if opt == 1,
   tz = griddata(X(:,1), X(:,2), X(:,3), tx, ty');
elseif opt == 2,
   tz = griddata(X(:,1), X(:,2), X(:,3), tx, ty','cubic');
elseif opt == 3,
   tz = griddata(X(:,1), X(:,2), X(:,3), tx, ty','nearest');
elseif opt == 4,
   tz = griddata(X(:,1), X(:,2), X(:,3), tx, ty','invdist');
end;
zmax = max(X(:,3));
zmin = min(X(:,3));
zdif = zmax - zmin;
lzd = floor(log10(zdif));
cdif = 10^(lzd-1);
s = 30;
cdif = cdif/2;
while s > 15
   cdif = 2*cdif;
   cmin = cdif*ceil(zmin/cdif);
   cmax = cdif*floor(zmax/cdif);
   v = [cmin:cdif:cmax];
   s = length(v);
end;
[C,h] = contour(tx,ty',tz, v);
clabel(C,h);
set(gca,'LineWidth',2);
hold off
figure
mesh(tx,ty,tz)
