%script to explore the grid and search radius problem
%jra, 9/2/06

xmin = 0;
ymin = 0;
xmax = 2;
ymax = 2;
%use grid_resolution = sqrt(2/pi); for search area = 1; it has a search
%readius of sqrt(1/pi) = 0.5642
grid_resolution = 1;

load r144.txt
% x = r958(:,1);
% y = r958(:,2);
% z = r958(:,3);
x = r144(:,1);
y = r144(:,2);
z = r144(:,3);
meanx = mean(x)
meany = mean(y)
xx=x-meanx;
yy=y-meany;

locs = find((xx >-2) & (xx<4) & (yy>-1.5) & (yy<3.5));


x = xmin:grid_resolution:xmax;
y = ymin:grid_resolution:ymax;
[X,Y]=meshgrid(x,y);

figure(1)
clf

subplot(2,2,1)
plot(X,Y, 'k+')
hold on
search_radius = grid_resolution.*sqrt(2);
for i = 1:length(x)
    for j = 1:length(y)
        circle(x(i),y(j),search_radius)
    end
end
axis equal
title('Search radius = grid resolution * sqrt(2)')
grid_area = grid_resolution.^2;
search_area = pi.*(search_radius.^2);
% gtext(['grid area = ',num2str(grid_area),' and search area = ',num2str(search_area), ' and search radius =', num2str(search_radius)]) 
text(-3, -2, ['grid area = ',num2str(grid_area),' and search area = ',num2str(search_area), ' and search radius =', num2str(search_radius)]) 
axis ([-2 4 -1.5 3.5])
xxx=xx(locs);
yyy=yy(locs);
plot(xxx,yyy,'k.')




subplot(2,2,2)
plot(X,Y, 'k+')
hold on
search_radius = grid_resolution;
for i = 1:length(x)
    for j = 1:length(y)
        circle(x(i),y(j),search_radius)
    end
end
axis equal
title('Search radius = grid resolution')
grid_area = grid_resolution.^2;
search_area = pi.*(search_radius.^2);
% gtext(['grid area = ',num2str(grid_area),' and search area = ',num2str(search_area), ' and search radius =', num2str(search_radius)]) 
text(-3, -2, ['grid area = ',num2str(grid_area),' and search area = ',num2str(search_area), ' and search radius =', num2str(search_radius)]) 
axis ([-2 4 -1.5 3.5])
plot(xxx,yyy,'k.')


subplot(2,2,4)
plot(X,Y, 'k+')
hold on

search_radius= sqrt((grid_resolution.^2)./pi);
for i = 1:length(x)
    for j = 1:length(y)
        circle(x(i),y(j),search_radius)
    end
end
axis equal
title('Search areas constant')
grid_area = grid_resolution.^2;
search_area = pi.*(search_radius.^2);
% gtext(['grid area = ',num2str(grid_area),' and search area = ',num2str(search_area), ' and search radius =', num2str(search_radius)]) 
text(-3, -2, ['grid area = ',num2str(grid_area),' and search area = ',num2str(search_area), ' and search radius =', num2str(search_radius)])  
axis ([-2 4 -1.5 3.5])
plot(xxx,yyy,'k.')


subplot(2,2,3)
plot(X,Y, 'k+')
hold on

search_radius= grid_resolution.*sqrt(2)./2;
for i = 1:length(x)
    for j = 1:length(y)
        circle(x(i),y(j),search_radius)
    end
end
axis equal
title('Search radius = grid resolution * sqrt(2) / 2')
grid_area = grid_resolution.^2;
search_area = pi.*(search_radius.^2);
% gtext(['grid area = ',num2str(grid_area),' and search area = ',num2str(search_area), ' and search radius =', num2str(search_radius)]) 
text(-3, -2, ['grid area = ',num2str(grid_area),' and search area = ',num2str(search_area), ' and search radius =', num2str(search_radius)]) 
plot(xxx,yyy,'k.')

locs = find((xx >1.5) & (xx<3) & (yy>-1) & (yy<0.5));
xxx=xx(locs);
yyy=yy(locs)
plot(xxx,yyy,'r.')
elev = z(locs);
for i=1:length(xxx)
     s = sprintf('%0.5g',elev(i));
     text(xxx(i)+0.1,yyy(i)+0.1,s)
end
axis ([-2 4 -1.5 3.5])