%sbd_cep_find.m

clear

data = [ 36.79140 -121.94833 3;
    36.77307 -121.92462 3;
    36.74418 -121.92462 4;
    36.77997 -121.98637 2;
    36.80884 -121.98637 3;
    36.78449 -121.88656 3;
    36.77997 -121.98637 3;
    36.78683 -122.04816 4;
    36.80283 -121.91026 2;
    36.77307 -121.92462 4;
    36.79140 -121.94833 3;
    36.77997 -121.98637 3;
    36.77997 -121.98637 6;
    36.78449 -121.88656 4;
    36.82028 -121.94833 3;
    36.79140 -121.94833 3;
    36.83171 -121.91026 4;
    ]

%plot(data(:,2),data(:,1),'o')

%find max/min and set up grid
lonmax = max(data(:,2));
lonmin = min(data(:,2));
latmax = max(data(:,1));
latmin = min(data(:,1));

gridlonmax = round(lonmax,5,'significant') + 0.02;
gridlonmin = round(lonmin,5,'significant') - 0.01;
gridlatmax = round(latmax,4,'significant') + 0.01;
gridlatmin = round(latmin,4,'significant') - 0.01;

lon = (gridlonmin:0.001:gridlonmax);
lat = (gridlatmin:0.001:gridlatmax);

heatgrid = zeros(size(lon,2),size(lat,2))';

for i = 1:size(lat,2)
    for j = 1:size(lon,2)
        for k = 1:size(data,1)
            heatgrid(i,j) = heatgrid(i,j) + isincircle( lon(j), lat(i), data(k,:));
        end
    end %end going through y
end %end going through x

s = surf(lon,lat,heatgrid);
s.EdgeColor = 'none';
view(2);
hold on
z = max(max(heatgrid)) +1;
for i=1:size(data,1)
    pos(i) = z;
end
plot3(data(:,2),data(:,1),pos,'ko');