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');
%       u{i}  = readvar('u.dat',XB(i),'2D');
%     u{i}  = u{i}/sqrt(sf(i));
% runup
    
    cd ..
end

% scale simulation time
% use random = 0
% scaled fneq, rt & dtbc   
% for scaled results output times should be equal
% use instat = 3 instead of instat = 4; we use PM spectrum  at depthscale 6
% scale taper time
% scaled hswitch (nd), dzmax (nd^1.5), hmin (nd), eps (nd)
xm = load('x.asc');
zm = load('z.asc');
time = load('t.asc');time=time.*3600;time(1)=1;
WL=0.806;%[4.2 0.806 0.585 0.461]
retreat=1.15; %[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
%    hline(ru10(i,2),':k'); 
   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}(1,2,1),'g-'); 
%       plot(x{1},u{1}(:,2,i)+zs{1}(1,2,1),'c-'); 
%    plot(x{2},u{2}(:,2,i)+zs{1}(1,2,1),'c--');
   title(num2str(i)); axis([87 99 0.4 1.3]);%axis([140 215 1 7]);%axis([86.6 88.6 0.25 0.7]);axis([-1 1 0.2 0.4]);
   pause(0.05); hold off;
end