clear all
close all
clc

grid_res = 5;
load('MAT\sb50m_data');
load('D:\werf_je\Werk\1.201.819-Proefsuppletie_Galgeplaat_2e_evaluatie\Data_analyse\Cirkel_bestek\cb_data');
titlenum = {'a','b','c','d'};

radius = (cb_data.opp/pi)^0.5;
xorg1 = max(cb_data.x)-radius;
xorg2 = min(cb_data.x)+radius;
xorg = mean([xorg1,xorg2]);
yorg1 = max(cb_data.y)-radius;
yorg2 = min(cb_data.y)+radius;
yorg = mean([yorg1,yorg2]);

figure
plot(cb_data.x,cb_data.y,'k-')
grid on
hold on
axis equal
H1=circle([xorg,yorg],radius+10,1000,':');
X10 = get(H1,'XData');
Y10 = get(H1,'YData');
%%
H2=circle([xorg,yorg],radius+50,1000,':');
X50 = get(H2,'XData');
Y50 = get(H2,'YData');
%%
H3=circle([xorg,yorg],radius+100,1000,':');
X100 = get(H3,'XData');
Y100 = get(H3,'YData');
%pdecirc(xc,yc,radius)

figure
for i = 1:length(sb50m_data.naam)
    clear X Y IN
    figure
    [C,h] = contourf(sb50m_data.x{i},sb50m_data.y{i},sb50m_data.z{i});
    hold on
    axis equal
    set(h,'LevelList',[-12:1:1])
    axis([0.530 0.570 3.960 4.000]*1e5)
    plot(cb_data.x,cb_data.y,'k-','LineWidth',1.0)
    xlabel('xRD (m)')
    ylabel('yRD (m)') 
    title(['(',titlenum{i},') bodemligging [m NAP], ',datestr(sb50m_data.datum{i}),', singlebeam 50 m'])
    caxis([-12 1])
    hcb = colorbar;
    set(hcb,'ytick',[-12:1:1])
    saveas(gcf,['EPS\sb50m_',num2str(i,'%02.0f'),'.eps'],'psc2')
    saveas(gcf,['PNG\sb50m_',num2str(i,'%02.0f'),'.png'])
    %Kubering
    [X,Y] = meshgrid(sb50m_data.x{i},sb50m_data.y{i});
    [IN] = inpolygon(X,Y,cb_data.x,cb_data.y);
    vol(i) = nansum(nansum(IN.*sb50m_data.z{i}))*grid_res^2;    
    %%%
    clear IN
    [IN] = inpolygon(X,Y,X10,Y10);
    vol_10(i) = nansum(nansum(IN.*sb50m_data.z{i}))*grid_res^2 - vol(i);        
    %%%
    clear IN
    [IN] = inpolygon(X,Y,X50,Y50);
    vol_50(i) = nansum(nansum(IN.*sb50m_data.z{i}))*grid_res^2;   
    %%%
    clear IN
    [IN] = inpolygon(X,Y,X100,Y100);
    vol_100(i) = nansum(nansum(IN.*sb50m_data.z{i}))*grid_res^2;       
end
sb50m_vol(:,1) = [sb50m_data.datum{:}];
sb50m_vol(:,2) = vol - vol(1);
sb50m_vol(:,3) = vol_10 - vol_10(1);
sb50m_vol(:,4) = vol_50 - vol_50(1);
sb50m_vol(:,5) = vol_100 - vol_100(1);

save MAT\sb50m_vol sb50m_vol