clear all
close all
warning off all
clc

load('MAT\rtkprof_data');
load('MAT\set_data');
load('MAT\sb25m_data');
load('MAT\rtk_data');
ioi = 1:2:7;
grp = {[90:93,101],[94:96],[97:100,101:102]};
ind1 ={[NaN,2,3,6,9,13,16],[NaN,2,3,6,9,13,16],[NaN,1,2,5,8,12,15],[NaN,2,3,6,9,13,16],[NaN,2,3,6,9,13,16],...
       [NaN,2,3,6,9,13,16],[NaN,2,3,6,9,13,16],[NaN,2,3,6,9,13,16],[NaN,2,3,6,9,13,16],[NaN,1,2,5,8,12,15],...
       [NaN,2,3,6,9,13,16],[NaN,NaN,NaN,1,4,8,11],[NaN,NaN,NaN,1,4,8,11],[NaN,NaN,NaN,1,4,8,11]};
ind2 = [NaN,3:4,8,NaN,NaN,NaN];
ind3 = [NaN,NaN,NaN,1,2,NaN,3];
dirr = {'(ZW - NO)','(W - O)','(ZO - NW)'};

for i = 1:length(rtkprof_data.naam)
    clear ind dx dy dist alfa beta gamma set_x X Y
    for j = 1:length(grp{i})
        clear dum1 dum2
        dum1 = [grp{i}];
        dum2 = strcmp(set_data.naam,num2str(dum1(j),'%03.0f'));
        ind(j) = find(dum2 == 1);
        dx(j) = set_data.x(ind(j)) - rtkprof_data.X0(i);
        dy(j) = set_data.y(ind(j)) - rtkprof_data.Y0(i);
        dist(j) = (dx(j)^2 + dy(j)^2)^0.5;             
    end
    beta = atan2(dy,dx)/pi*180;
    alfa = rtkprof_data.theta(i);
    gamma = (270 - alfa)- beta;
    set_x = -cos(gamma/180*pi).*dist;    
    for k = 1:length(rtkprof_data.datum{i})
        clear dum1 dum2 dum4 sb25m_z sb25m_z2 xRD yRD rtk_z rtk_z2 X2 Y2
        dum1 = rtkprof_data.X{i};
        dum2 = rtkprof_data.Z{i}; 
        xRD = rtkprof_data.xRD{i};
        yRD = rtkprof_data.yRD{i};
        if ~isnan(ind2(k))
            sb25m_z = sb25m_data.z{ind2(k)};
            [X,Y] = meshgrid(sb25m_data.x{ind2(k)},sb25m_data.y{ind2(k)});
            sb25m_z2 = interp2(X, Y, sb25m_z, xRD{k}, yRD{k});
        else
            X = NaN;
            Y = NaN;
            sb25m_z = NaN*dum1{k};            
            sb25m_z2 = NaN*dum1{k};
        end
        if ~isnan(ind3(k))
            rtk_z = rtk_data.z{ind3(k)};
            [X2,Y2] = meshgrid(rtk_data.x{ind3(k)},rtk_data.y{ind3(k)});
            rtk_z2 = interp2(X2, Y2, rtk_z, xRD{k}, yRD{k});
        else
            X2 = NaN;
            Y2 = NaN;
            rtk_z = NaN*dum1{k};            
            rtk_z2 = NaN*dum1{k};
        end
            
        for j = 1:length(grp{i})
            clear dum3
            dum3 = set_data.z{ind(j)};
            id1 = ind1{ind(j)};
            if ~isnan(id1(k))
                dum4(:,j) = [dum3(id1(k))];
            else
                dum4(:,j) = NaN;
            end
        end
        figure
        subplot(2,1,1)
        plot(dum1{k},dum2{k},'k-','LineWidth',1.0)
        hold on
        plot(set_x,dum4,'r^','MarkerFaceColor','r')
        plot(dum1{k},sb25m_z2,'b-','LineWidth',1.0)
        plot(dum1{k},rtk_z2,'g-','LineWidth',1.0)
        grid on                
        xlim([-400 400])     
        plot(xlim,xlim.*0,'k-')   
        dumm = rtkprof_data.datum{i};
        title([rtkprof_data.naam{i},' ',dirr{i},', ',datestr(dumm(k))])        
        ylabel('z_b (m NAP)')
        if i == 3
            legend('RTK profiel','SET','singlebeam 25 m','RTK',2)
            xlabel('positie langs profiel (m)')
        else
            legend('RTK profiel','SET','singlebeam 25 m','RTK',1)
        end        
        saveas(gcf,['PNG\comp-rtk_prof_',num2str(i,'%02.0f'),'_',num2str(k,'%02.0f'),'.png'])
        saveas(gcf,['EPS\comp-rtk_prof_',num2str(i,'%02.0f'),'_',num2str(k,'%02.0f'),'.eps'],'psc2')
    end  
end
