clear all
close all
clc

load('dvd_data.mat');

grid_res = 2.5;
load('D:\werf_je\Werk\1.201.819-Proefsuppletie_Galgeplaat_2e_evaluatie\Data_analyse\Bodemligging\rtk_data');
load('D:\werf_je\Werk\1.201.819-Proefsuppletie_Galgeplaat_2e_evaluatie\Data_analyse\Bodemligging\sb25m_data');
load('D:\werf_je\Werk\1.201.819-Proefsuppletie_Galgeplaat_2e_evaluatie\Data_analyse\Cirkel_bestek\cb_data');

%Gecombineerde set maken
len1 = length(rtk_data.datum);
len2 = length(sb25m_data.datum);
%
rtk_sb25m_data.datum(1:len1) = rtk_data.datum;
rtk_sb25m_data.datum(len1+1:len1+len2) = sb25m_data.datum;
%
rtk_sb25m_data.naam(1:len1) = rtk_data.naam;
rtk_sb25m_data.naam(len1+1:len1+len2) = sb25m_data.naam;
%
rtk_sb25m_data.x(1:len1) = rtk_data.x;
rtk_sb25m_data.x(len1+1:len1+len2) = sb25m_data.x;
%
rtk_sb25m_data.y(1:len1) = rtk_data.y;
rtk_sb25m_data.y(len1+1:len1+len2) = sb25m_data.y;
%
rtk_sb25m_data.z(1:len1) = rtk_data.z;
rtk_sb25m_data.z(len1+1:len1+len2) = sb25m_data.z;
%
rtk_sb25m_data.type(1:len1) = {'RTK'};
rtk_sb25m_data.type(len1+1) = {'singlebeam 25 m'};
rtk_sb25m_data.type(len1+2:len1+4) = {'gecombineerde meting'};
rtk_sb25m_data.type(len1+5:len1+len2) = {'singlebeam 25 m'};

%Even sorten op datum
[hm ind] = sort([rtk_sb25m_data.datum{:}]);
rtk_sb25m_data.datum = rtk_sb25m_data.datum(ind);
rtk_sb25m_data.naam = rtk_sb25m_data.naam(ind);
rtk_sb25m_data.x = rtk_sb25m_data.x(ind);
rtk_sb25m_data.y = rtk_sb25m_data.y(ind);
rtk_sb25m_data.z = rtk_sb25m_data.z(ind);
rtk_sb25m_data.type= rtk_sb25m_data.type(ind);

levels = dvd_data.levels{:};

h = waitbar(0,'Metingen...');
for i = 1:length(rtk_sb25m_data.naam)
    waitbar(i/length(rtk_sb25m_data.naam),h);
    clear X Y Z
    [X,Y] = meshgrid(rtk_sb25m_data.x{i},rtk_sb25m_data.y{i});
    Z = rtk_sb25m_data.z{i};
    [IN] = inpolygon(X,Y,cb_data.x,cb_data.y);
    h2 = waitbar(0,'Niveaus...');
    for j = 1:length(levels)
        waitbar(j/length(levels),h2);
        clear dum
        dum = find(Z(IN) > levels(j));
        if ~isempty(dum)
            areaal(j,i) = length(dum)*grid_res^2;
        else
            areaal(j,i) = 0;
        end            
    end
    close(h2);
    %Maken data struct
    dva_data.areaal(i) = {areaal(:,i)};
    dva_data.datum(i) = rtk_sb25m_data.datum(i);
    dva_data.naam(i) = rtk_sb25m_data.naam(i);    
    if i == 5 | i == 9 | i == 13
        figure        
        plot(levels,dva_data.areaal{1},'k-','LineWidth',1)
        hold on
        grid on
        plot(levels,dva_data.areaal{i-3},'r-','LineWidth',1)
        plot(levels,dva_data.areaal{i-2},'b-','LineWidth',1)
        plot(levels,dva_data.areaal{i-1},'g-','LineWidth',1)
        plot(levels,dva_data.areaal{i},'m-','LineWidth',1)
        xlabel('z_w (m NAP)')
        ylabel('droogvallend areaal (m^2)')
        legend(datestr(dva_data.datum{1}),datestr(dva_data.datum{i-3}),datestr(dva_data.datum{i-2}),...
               datestr(dva_data.datum{i-1}),datestr(dva_data.datum{i})); 
        xlim([-1.5 0.5])
        pause        
    end
end
close(h);

%Wegschrijven data
save dva_data dva_data