close all; clear all;

addpath P:\1002266-xbeach\code

dirs = {'run000'};
eval(['cd ',dirs{1}]); 
XB(1) = getdimensions;
cd ..

for i = 1:length(XB)
    
    x{i}  = XB(i).xc;
    
    eval(['cd ',dirs{i}]); 
    zb{i} = readvar('zb.dat',XB(i),'2D');    
    zs{i} = readvar('zs.dat',XB(i),'2D');  
    H{i}  = readvar('H.dat',XB(i),'2D');
    cd ..
end

xm = load('x.asc');
zm = load('z.asc');
time = load('t.asc');time=time.*3600;time(1)=1;
WL=0.461;%[4.2 0.806 0.585 0.461]
retreat=0.58; %[5.8 1.15 0.78 0.58]
for i=1:length(time); 
    BSS(i) = BrierSkillScore(x{1}(:,2), zb{1}(:,2,time(i)), xm(:,i), zm(:,i), xm(:,1), zm(:,1),10000);
    [Volume result Boundaries] = getVolume('x2',x{1}(:,2), 'z2', zb{1}(:,2,time(i)), 'x',x{1}(:,2),'z',zb{1}(:,2,1),'LowerBoundary',WL);
    volumeXB(i)=result.Volumes.Erosion;
    [Volume result Boundaries] = getVolume('x2',xm(:,i), 'z2', zm(:,i), 'x',xm(:,1),'z',zm(:,1),'LowerBoundary',WL);
    volumeM(i)=result.Volumes.Erosion;
    left  = findCrossings(x{1}(:,2),zb{1}(:,2,1),x{1}(:,2),repmat(retreat,1,length(x{1}(:,2)))');
    right = findCrossings(x{1}(:,2),zb{1}(:,2,time(i)),x{1}(:,2),repmat(retreat,1,length(x{1}(:,2)))');
    retreatXB(i) = abs(left - right);
    left  = findCrossings(xm(~isnan(xm(:,1)),1),zm(~isnan(xm(:,1)),1),xm(~isnan(xm(:,1)),1),repmat(retreat,1,length(xm(~isnan(xm(:,1)),1)))');
    right = findCrossings(xm(~isnan(xm(:,i)),i),zm(~isnan(xm(:,i)),i),xm(~isnan(xm(:,i)),i),repmat(retreat,1,length(xm(~isnan(xm(:,i)),i)))');
    retreatM(i) = abs(left - right);
end
BSS,volumeXB,volumeM,retreatXB,retreatM

figure; 
set(gcf,'position',[1 1 1200 860])
for i = 1:100:XB(1).nt
   plot(xm,zm,':k'),hold on
   plot(x{1},zb{1}(:,2,1),'k--'); 
   plot(x{1},zb{1}(:,2,i),'k-o'); 
   plot(x{1},zs{1}(:,2,i),'b-'); 
   plot(x{1},H{1}(:,2,i)+zs{1}(:,2,i),'g-'); 
   title(num2str(i)); axis([73 79 0.25 0.7]);%axis([92.5 96 0.3 0.9]);%axis([140 215 1 7]);%axis([86.6 88.6 0.25 0.7]);axis([-1 1 0.2 0.4]);
   legend('Measured','Initial','Bed level','Water level','Wave height',2)
   pause(0.05); hold off;
end

j=3;
figure(),hold on,box on
        plot(xm(:,1),zm(:,1),'-ok','MarkerSize',3);
        plot(xm(:,j),zm(:,j),':ok','MarkerSize',3);
        plot(x{1},zb{1}(:,2,21600),':xr','MarkerSize',5);
        legend('M 0.0h','M 6.0h','XBeach 6.0h',4);
        xlabel('X-direction from waveboard [m]'), ylabel('Elevation [m]')
        axis([84.5 88.8 0.31 0.61])
        hline(0.461,':b','WL')