%% Simple example show the Volumetric changes in wet area
clear
%close all
clc

savefigs=1;

%% control zones

% Cross shore
B_offs_all= 1200;
B_ons_all = 0;
B_offs_1  = 100;
B_ons_1   = 0;
B_offs_2  = 300;
B_ons_2   = 100;
B_offs_3  = 500;
B_ons_3   = 300;
B_offs_4  = 800 ;
B_ons_4   = 500;
B_offs_5  = 1200;
B_ons_5   = 800;

%%  load survey data

Mpath='F:\check_outs\tudelft\vlugtenburg\';
%Mpath='C:\Madeschipper\check_out\OpenEarthRawData\tudelft\vlugtenburg\scripts\Matthieu'

load([Mpath '\AllVbeachsurveys.mat'])
load([Mpath '\datesandpairs.mat'])

disp('data loaded')

%%
%% Fix data
disp('--- precheck')
% find spots where there is no data at or befor the x=0 point
for is=1:length(Vsurvey);
    disp(['survey  ' num2str(is)])
    for it=1:22
        if isnan(Vsurvey(is,1).transects.elevation2(301:311,it))
            disp(['error.. no data near -9 m NAP ' num2str(it)])
        end
        if isnan(Vsurvey(is,1).transects.elevation2(51:61,it))
            disp(['error.. no data near dune fence ' num2str(it)])
        end
        
    end
end

% Fill data with surveys from before/after
%%
Vsurvey(8,1).transects.elevation2(61,21)=Vsurvey(7,1).transects.elevation2(61,21);
Vsurvey(8,1).transects.elevation2(61,22)=Vsurvey(7,1).transects.elevation2(61,22);
%%
Vsurvey(1,1).transects.elevation2(301:311,1:2)=Vsurvey(2,1).transects.elevation2(301:311,1:2);
Vsurvey(1,1).transects.elevation2(301:311,1:2)=Vsurvey(2,1).transects.elevation2(301:311,1:2);
Vsurvey(4,1).transects.elevation2(301:311,1:4)=Vsurvey(3,1).transects.elevation2(301:311,1:4);
Vsurvey(4,1).transects.elevation2(301:311,1:4)=Vsurvey(3,1).transects.elevation2(301:311,1:4);
Vsurvey(9,1).transects.elevation2(51:61,22)=Vsurvey(8,1).transects.elevation2(51:61,22);
Vsurvey(10,1).transects.elevation2(301:311,1)=Vsurvey(9,1).transects.elevation2(301:311,1);
Vsurvey(11,1).transects.elevation2(51:61,22)=Vsurvey(10,1).transects.elevation2(51:61,22);
Vsurvey(12,1).transects.elevation2(51:61,22)=Vsurvey(13,1).transects.elevation2(51:61,22);
Vsurvey(14,1).transects.elevation2(51:61,22)=Vsurvey(13,1).transects.elevation2(51:61,22);
Vsurvey(15,1).transects.elevation2(301:311,1:2)=Vsurvey(14,1).transects.elevation2(301:311,1:2);
Vsurvey(19,1).transects.elevation2(301:311,1:2)=Vsurvey(20,1).transects.elevation2(301:311,1:2);
Vsurvey(18,1).transects.elevation2(301:311,1:8)=Vsurvey(19,1).transects.elevation2(301:311,1:8);
Vsurvey(17,1).transects.elevation2(301:311,1:5)=Vsurvey(18,1).transects.elevation2(301:311,1:5);
Vsurvey(22,1).transects.elevation2(51:61,22)=Vsurvey(21,1).transects.elevation2(51:61,22);
Vsurvey(25,1).transects.elevation2(51:61,22)=Vsurvey(24,1).transects.elevation2(51:61,22);
Vsurvey(28,1).transects.elevation2(51:61,22)=Vsurvey(27,1).transects.elevation2(51:61,22);
Vsurvey(31,1).transects.elevation2(51:61,22)=Vsurvey(30,1).transects.elevation2(51:61,22);

Vsurvey(35,1).transects.elevation2(51:61,22)=Vsurvey(34,1).transects.elevation2(51:61,22);

for is=1:length(Vsurvey);
    disp(['survey  ' num2str(is)])
    for it=1:22
        if isnan(Vsurvey(is,1).transects.elevation2(301:311,it))
            disp(['error.. no data near -9 m NAP ' num2str(it)])
        end
        if isnan(Vsurvey(is,1).transects.elevation2(51:61,it))
            disp(['error.. no data near dune fence ' num2str(it)])
        end
    end
end



%%

up=1.5;
down=-4;

dz=0.1;

b_b=up-0.5*dz;
b_up=up+0.5*dz;

b_b2=down-0.5*dz;
b_up2=down+0.5*dz;

B_offs_all= 1200;
B_ons_all = 0;

  for i_date=1:length(Vsurvey)
      for i_transect=1:22;
            
            xi=Vsurvey(i_date).transects.cr_sh;
            zi=Vsurvey(i_date).transects.elevation2(:,i_transect);
            Vol_up(i_date,i_transect)=jarkus_getVolume(xi, zi, b_up, b_b ,B_ons_all,B_offs_all);
     x_up(i_date,i_transect)=Vol_up(i_date,i_transect)/(b_up-b_b);
     
                 Vol_down(i_date,i_transect)=jarkus_getVolume(xi, zi, b_up2, b_b2 ,B_ons_all,B_offs_all);
     x_down(i_date,i_transect)=Vol_down(i_date,i_transect)/(b_up2-b_b2);
     
     
     
     if i_date<=23 && i_transect==11
        x_up(i_date,i_transect)=NaN;
     end
       slope(i_date,i_transect)=(x_down(i_date,i_transect)-x_up(i_date,i_transect))/(up - down);
      end
  end

disp('slope calculated!')
  
  
  
%%  fig prep

monthv=[];
for i_mm=7:12
    iv=datenum([num2str(i_mm) '/01/2009'],'mm/dd/yyyy');
    monthv=[ monthv iv];
end
for i_mm=1:12
    iv=datenum([num2str(i_mm) '/01/2010'],'mm/dd/yyyy');
    monthv=[ monthv iv];
end
for i_mm=1:12
    iv=datenum([num2str(i_mm) '/01/2011'],'mm/dd/yyyy');
    monthv=[ monthv iv];
end
for i_mm=1:11
    iv=datenum([num2str(i_mm) '/01/2012'],'mm/dd/yyyy');
    monthv=[ monthv iv];
end

xtickstring=datestr(monthv,'mmm/yy');
xtickstring(1:2:end,:)=NaN;

colors=colormap(winter(22));

%%
figure
set(gcf,'position',[  50         203        1310         378])
set(gcf,'PaperPosition',[0 0  35 12])

hold on
for i_transect=1:22
    for i_date=1:1:length(Vsurvey)
        datesmatlab(i_date)=Vsurvey(i_date).matlabdate;
plot(datesmatlab(i_date),slope(i_date,i_transect),'.','markersize',14,'color',colors(i_transect,:))
% plot(i_date,Bar(i_date,1).barposition2(i_transect,1),'.r')

    end
end
colormap(colors);
h=colorbar;
     
h2=annotation('textarrow',[0.92 0.7],[0.13 0.8],'string',{'Transect nr. , 1 (North) to 22 (South)'}, ...
            'HeadStyle','none','LineStyle', 'none', 'TextRotation',-90);
% 
    for i_date=1:1:length(Vsurvey)
%         barsurvtemp=[]
%         for i_transect=1:22
%             barsurvtemp(i_transect)=Bar(i_date,1).barposition(i_transect,1)
% % plot(i_date,Bar(i_date,1).barposition2(i_transect,1),'.r')
% 
%         end
     MeanSlope(i_date)=mean(slope(i_date,~isnan(slope(i_date,:))))  
     end
% 
 plot(datesmatlab,MeanSlope,'color',[0.5 0 0])

xlim([datenum(['06/01/2009'],'mm/dd/yyyy') datenum(['11/1/2012'],'mm/dd/yyyy')])

xlabel('Date')
ylabel({'Cross-shore slope';' between -2 and -4 m NAP [1: ]'})


holdax=axis;

x = [datenum(['01/01/2010'],'mm/dd/yyyy'); datenum(['12/31/2010'],'mm/dd/yyyy');datenum(['12/31/2010'],'mm/dd/yyyy') ; datenum(['01/01/2010'],'mm/dd/yyyy')];
y = [ 600; 600;0; 0];
tcolor(1,1,1:3) = 0.7*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')

x = [datenum(['01/01/2012'],'mm/dd/yyyy'); datenum(['12/31/2012'],'mm/dd/yyyy');datenum(['12/31/2012'],'mm/dd/yyyy') ; datenum(['01/01/2012'],'mm/dd/yyyy')];
y = [ 600; 600;0; 0];
tcolor(1,1,1:3) = 0.7*[1 1 1];
patch(x,y,tcolor,'facealpha',0.3,'EdgeColor','none')


text(datenum(['06/01/2010'],'mm/dd/yyyy'),15,'2010','color','w','fontsize',18,'FontWeight','bold')
text(datenum(['04/10/2012'],'mm/dd/yyyy'),15,'2012','color','w','fontsize',18,'FontWeight','bold')
text(datenum(['09/01/2009'],'mm/dd/yyyy'),15,'2009','color',0.7*[1 1 1],'fontsize',18,'FontWeight','bold')
text(datenum(['06/01/2011'],'mm/dd/yyyy'),15,'2011','color',0.7*[1 1 1],'fontsize',18,'FontWeight','bold')
grid on;

axis(holdax);
set(gca,'Xtick',monthv)
set(gca,'Xticklabel',xtickstring)

ylim([0 140])

set (gcf,'color','w')
%drawnow
print ([ 'ProfileSlope.png'],'-dpng')
