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

savefigs=0;

%% control zones

% Cross shore
B_offs_all = 1200;
B_ons_all = 000;
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\scripts\Matthieu';
%Mpath='C:\Madeschipper\check_out\OpenEarthRawData\tudelft\vlugtenburg\scripts\Matthieu'

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

disp('data loaded')

%% Fix
Vsurvey(8,1).transects.elevation(61,21)=Vsurvey(7,1).transects.elevation(61,21);
Vsurvey(8,1).transects.elevation(61,22)=Vsurvey(7,1).transects.elevation(61,22);
%
Vsurvey(1,1).transects.elevation(301:311,1:2)=Vsurvey(2,1).transects.elevation(301:311,1:2);
Vsurvey(1,1).transects.elevation(301:311,1:2)=Vsurvey(2,1).transects.elevation(301:311,1:2);
Vsurvey(4,1).transects.elevation(301:311,1:4)=Vsurvey(3,1).transects.elevation(301:311,1:4);
Vsurvey(4,1).transects.elevation(301:311,1:4)=Vsurvey(3,1).transects.elevation(301:311,1:4);
Vsurvey(9,1).transects.elevation(51:61,22)=Vsurvey(8,1).transects.elevation(51:61,22);
Vsurvey(10,1).transects.elevation(301:311,1)=Vsurvey(9,1).transects.elevation(301:311,1);
Vsurvey(11,1).transects.elevation(51:61,22)=Vsurvey(10,1).transects.elevation(51:61,22);
Vsurvey(12,1).transects.elevation(51:61,22)=Vsurvey(13,1).transects.elevation(51:61,22);
Vsurvey(14,1).transects.elevation(51:61,22)=Vsurvey(13,1).transects.elevation(51:61,22);
Vsurvey(15,1).transects.elevation(301:311,1:2)=Vsurvey(14,1).transects.elevation(301:311,1:2);
Vsurvey(19,1).transects.elevation(301:311,1:2)=Vsurvey(20,1).transects.elevation(301:311,1:2);
Vsurvey(18,1).transects.elevation(301:311,1:8)=Vsurvey(19,1).transects.elevation(301:311,1:8);
Vsurvey(17,1).transects.elevation(301:311,1:5)=Vsurvey(18,1).transects.elevation(301:311,1:5);
Vsurvey(22,1).transects.elevation(51:61,22)=Vsurvey(21,1).transects.elevation(51:61,22);
Vsurvey(25,1).transects.elevation(51:61,22)=Vsurvey(24,1).transects.elevation(51:61,22);
Vsurvey(28,1).transects.elevation(51:61,22)=Vsurvey(27,1).transects.elevation(51:61,22);
%
% fix SoS
Vsurvey(7,1).transects.elevation(Vsurvey(7,1).transects.elevation <0)=Vsurvey(7,1).transects.elevation(Vsurvey(7,1).transects.elevation <0)*0.990;
Vsurvey(10,1).transects.elevation(Vsurvey(10,1).transects.elevation <0)=Vsurvey(10,1).transects.elevation(Vsurvey(10,1).transects.elevation <0)*1.005;
Vsurvey(12,1).transects.elevation(Vsurvey(12,1).transects.elevation <0)=Vsurvey(12,1).transects.elevation(Vsurvey(12,1).transects.elevation <0)*0.995;
Vsurvey(13,1).transects.elevation(Vsurvey(13,1).transects.elevation <0)=Vsurvey(13,1).transects.elevation(Vsurvey(13,1).transects.elevation <0)*1.005;
Vsurvey(14,1).transects.elevation(Vsurvey(14,1).transects.elevation <0)=Vsurvey(14,1).transects.elevation(Vsurvey(14,1).transects.elevation <0)*0.995;
Vsurvey(16,1).transects.elevation(Vsurvey(16,1).transects.elevation <0)=Vsurvey(16,1).transects.elevation(Vsurvey(16,1).transects.elevation <0)*1.005;
Vsurvey(17,1).transects.elevation(Vsurvey(17,1).transects.elevation <0)=Vsurvey(17,1).transects.elevation(Vsurvey(17,1).transects.elevation <0)*0.995;
Vsurvey(18,1).transects.elevation(Vsurvey(18,1).transects.elevation <0)=Vsurvey(18,1).transects.elevation(Vsurvey(18,1).transects.elevation <0)*0.995;
Vsurvey(19,1).transects.elevation(Vsurvey(19,1).transects.elevation <0)=Vsurvey(19,1).transects.elevation(Vsurvey(19,1).transects.elevation <0)*0.995;




%% Jarkus method
% V vertical (planview boxes)
b_b=-10;
b_up=10;

dy=-1*diff(Vsurvey(1,1).transects.y);
for i_date=1:length(Vsurvey)
    disp(['' num2str(i_date) ])
    for i_transect=1:22;
        xi=Vsurvey(i_date).transects.cr_sh;
        zi=Vsurvey(i_date).transects.elevation(:,i_transect);
        Vsurvey(i_date).VolV(i_transect)=jarkus_getVolume(xi, zi, b_up, b_b ,B_ons_all,B_offs_all);
        
        Vsurvey(i_date).VolV1(i_transect)=jarkus_getVolume(xi, zi, b_up, b_b ,B_ons_1,B_offs_1);
        Vsurvey(i_date).VolV2(i_transect)=jarkus_getVolume(xi, zi, b_up, b_b ,B_ons_2,B_offs_2);
        Vsurvey(i_date).VolV3(i_transect)=jarkus_getVolume(xi, zi, b_up, b_b ,B_ons_3,B_offs_3);
        Vsurvey(i_date).VolV4(i_transect)=jarkus_getVolume(xi, zi, b_up, b_b ,B_ons_4,B_offs_4);
        Vsurvey(i_date).VolV5(i_transect)=jarkus_getVolume(xi, zi, b_up, b_b ,B_ons_5,B_offs_5);
       
        if i_transect==1
            dy_tr=0.5*dy(1);
        end
        if i_transect==22
            dy_tr=0.5*dy(21);
        end
        if i_transect>=2 && i_transect<=21
            dy_tr=0.5*dy(i_transect-1)+0.5*dy(i_transect);
        end
        
        Vsurvey(i_date).transects.dy_tr(i_transect)=dy_tr;
        
        Vsurvey(i_date).VolV_tr(i_transect)=Vsurvey(i_date).VolV(i_transect)*dy_tr;
        Vsurvey(i_date).VolV1_tr(i_transect)=Vsurvey(i_date).VolV1(i_transect)*dy_tr;
        Vsurvey(i_date).VolV2_tr(i_transect)=Vsurvey(i_date).VolV2(i_transect)*dy_tr;
        Vsurvey(i_date).VolV3_tr(i_transect)=Vsurvey(i_date).VolV3(i_transect)*dy_tr;
        Vsurvey(i_date).VolV4_tr(i_transect)=Vsurvey(i_date).VolV4(i_transect)*dy_tr;
        Vsurvey(i_date).VolV5_tr(i_transect)=Vsurvey(i_date).VolV5(i_transect)*dy_tr;
        
    end
    Vsurvey(i_date).Vtotal=sum(Vsurvey(i_date).VolV_tr);
    Vsurvey(i_date).Vtotal1=sum(Vsurvey(i_date).VolV1_tr);
    Vsurvey(i_date).Vtotal2=sum(Vsurvey(i_date).VolV2_tr);
    Vsurvey(i_date).Vtotal3=sum(Vsurvey(i_date).VolV3_tr);
    Vsurvey(i_date).Vtotal4=sum(Vsurvey(i_date).VolV4_tr);
    Vsurvey(i_date).Vtotal5=sum(Vsurvey(i_date).VolV5_tr);
    Vtotal(i_date)=sum(Vsurvey(i_date).VolV_tr);
    Vtotal1(i_date)=sum(Vsurvey(i_date).VolV1_tr);
    Vtotal2(i_date)=sum(Vsurvey(i_date).VolV2_tr);
    Vtotal3(i_date)=sum(Vsurvey(i_date).VolV3_tr);
    Vtotal4(i_date)=sum(Vsurvey(i_date).VolV4_tr);
    Vtotal5(i_date)=sum(Vsurvey(i_date).VolV5_tr);
end


%% Check is volumes add up
for i_date=1:length(Vsurvey)
    
    difference=Vtotal(i_date) - ( Vtotal1(i_date)+Vtotal2(i_date)+Vtotal3(i_date)+Vtotal4(i_date)+Vtotal5(i_date));
    if difference >=1
        i_date
        disp('Volumes do not add up!')
    end
end

%% Figure Volume changes wrt to first survey in boxes
% Planfom boxes!
figure

plot(datenum(dates(pairs(1,1),:)),0,'or') % zero point
hold on;
bar(datenum(dates(pairs(:,1),:)),(Vtotal-Vtotal(1))/1745)
plot(datenum(dates(pairs(:,1),:)),(Vtotal1-Vtotal1(1))/1745,'*--g')
plot(datenum(dates(pairs(:,1),:)),(Vtotal2-Vtotal2(1))/1745,'o--m')
plot(datenum(dates(pairs(:,1),:)),(Vtotal3-Vtotal3(1))/1745,'+--k')
plot(datenum(dates(pairs(:,1),:)),(Vtotal4-Vtotal4(1))/1745,'v--b')
plot(datenum(dates(pairs(:,1),:)),(Vtotal5-Vtotal5(1))/1745,'^--r')
%ylim([-250 200])
xlim([ datenum(dates(pairs(1,1),:)) datenum(dates(pairs(end,1),:)) ])
datetick('x')
grid on
ylabel('\Delta V per m alongshore wrt June 2009')
legend('First survey','Total area','Area I','Area II','Area III','Area IV','Area V','location','northwest')

set(gcf,'PaperPosition',[0  0  25 20])
hold off

title('Volume method, Boxes')

if savefigs==1
    filename=[ 'Volume control area y=  ' num2str(B_south) '_' num2str(B_north) ', x=  ' num2str(B_ons_all) '_' num2str(B_offs_all)]
    print ([ Mpath '\VolumeAnalysis\images\' filename '.jpg'],'-djpeg')
end


%% Jarkus method
% V horizontal (discs )

b_b=-10;
b_up=10;

b_bA=2.5;
b_upA=10;

b_bB =-1.5;
b_upB=2.5;

b_bC=-5;
b_upC=-1.5;

b_bD=-7;
b_upD=-5;

b_bE=-10;
b_upE=-7;


for i_date=1:length(Vsurvey)
    disp(['' num2str(i_date) ])
    for i_transect=1:22;
        xi=Vsurvey(i_date).transects.cr_sh;
        zi=Vsurvey(i_date).transects.elevation(:,i_transect);
        
        Vsurvey(i_date).VolVA(i_transect)=jarkus_getVolume(xi, zi, b_upA , b_bA ,B_ons_all,B_offs_all);
        Vsurvey(i_date).VolVB(i_transect)=jarkus_getVolume(xi, zi, b_upB, b_bB ,B_ons_all,B_offs_all);
        Vsurvey(i_date).VolVC(i_transect)=jarkus_getVolume(xi, zi, b_upC, b_bC ,B_ons_all,B_offs_all);
        Vsurvey(i_date).VolVD(i_transect)=jarkus_getVolume(xi, zi, b_upD, b_bD ,B_ons_all,B_offs_all);
        Vsurvey(i_date).VolVE(i_transect)=jarkus_getVolume(xi, zi, b_upE, b_bE  ,B_ons_all,B_offs_all);
       
        if i_transect==1
            dy_tr=0.5*dy(1);
        end
        if i_transect==22
            dy_tr=0.5*dy(21);
        end
        if i_transect>=2 && i_transect<=21
            dy_tr=0.5*dy(i_transect-1)+0.5*dy(i_transect);
        end
        
        Vsurvey(i_date).transects.dy_tr(i_transect)=dy_tr;
        
        Vsurvey(i_date).VolVA_tr(i_transect)=Vsurvey(i_date).VolVA(i_transect)*dy_tr;
        Vsurvey(i_date).VolVB_tr(i_transect)=Vsurvey(i_date).VolVB(i_transect)*dy_tr;
        Vsurvey(i_date).VolVC_tr(i_transect)=Vsurvey(i_date).VolVC(i_transect)*dy_tr;
        Vsurvey(i_date).VolVD_tr(i_transect)=Vsurvey(i_date).VolVD(i_transect)*dy_tr;
        Vsurvey(i_date).VolVE_tr(i_transect)=Vsurvey(i_date).VolVE(i_transect)*dy_tr;
        
    end

    Vsurvey(i_date).VtotalA=sum(Vsurvey(i_date).VolVA_tr);
    Vsurvey(i_date).VtotalB=sum(Vsurvey(i_date).VolVB_tr);
    Vsurvey(i_date).VtotalC=sum(Vsurvey(i_date).VolVC_tr);
    Vsurvey(i_date).VtotalD=sum(Vsurvey(i_date).VolVD_tr);
    Vsurvey(i_date).VtotalE=sum(Vsurvey(i_date).VolVE_tr);

    VtotalA(i_date)=sum(Vsurvey(i_date).VolVA_tr);
    VtotalB(i_date)=sum(Vsurvey(i_date).VolVB_tr);
    VtotalC(i_date)=sum(Vsurvey(i_date).VolVC_tr);
    VtotalD(i_date)=sum(Vsurvey(i_date).VolVD_tr);
    VtotalE(i_date)=sum(Vsurvey(i_date).VolVE_tr);
end
%% Check is volumes add up
for i_date=1:length(Vsurvey)
    
    difference=Vtotal(i_date) - ( VtotalA(i_date)+VtotalB(i_date)+VtotalC(i_date)+VtotalD(i_date)+VtotalE(i_date));
    if difference >=1
        i_date
        disp('Volumes do not add up!')
    end
end

%% Jarkus method extra deep!!
% V horizontal (discs )

b_b=-10;
b_up=10;

b_bA=-6.5;
b_upA=-6;

b_bB =-7;
b_upB=-6.5;

b_bC=-7.5;
b_upC=-7;

b_bD=-8;
b_upD=-7.5;

b_bE=-10;
b_upE=-8;


for i_date=1:length(Vsurvey)
    disp(['' num2str(i_date) ])
    for i_transect=1:22;
        xi=Vsurvey(i_date).transects.cr_sh;
        zi=Vsurvey(i_date).transects.elevation(:,i_transect);
        
        Vsurvey(i_date).VolVA(i_transect)=jarkus_getVolume(xi, zi, b_upA , b_bA ,B_ons_all,B_offs_all);
        Vsurvey(i_date).VolVB(i_transect)=jarkus_getVolume(xi, zi, b_upB, b_bB ,B_ons_all,B_offs_all);
        Vsurvey(i_date).VolVC(i_transect)=jarkus_getVolume(xi, zi, b_upC, b_bC ,B_ons_all,B_offs_all);
        Vsurvey(i_date).VolVD(i_transect)=jarkus_getVolume(xi, zi, b_upD, b_bD ,B_ons_all,B_offs_all);
        Vsurvey(i_date).VolVE(i_transect)=jarkus_getVolume(xi, zi, b_upE, b_bE  ,B_ons_all,B_offs_all);
       
        if i_transect==1
            dy_tr=0.5*dy(1);
        end
        if i_transect==22
            dy_tr=0.5*dy(21);
        end
        if i_transect>=2 && i_transect<=21
            dy_tr=0.5*dy(i_transect-1)+0.5*dy(i_transect);
        end
        
        Vsurvey(i_date).transects.dy_tr(i_transect)=dy_tr;
        
        Vsurvey(i_date).VolVA_tr(i_transect)=Vsurvey(i_date).VolVA(i_transect)*dy_tr;
        Vsurvey(i_date).VolVB_tr(i_transect)=Vsurvey(i_date).VolVB(i_transect)*dy_tr;
        Vsurvey(i_date).VolVC_tr(i_transect)=Vsurvey(i_date).VolVC(i_transect)*dy_tr;
        Vsurvey(i_date).VolVD_tr(i_transect)=Vsurvey(i_date).VolVD(i_transect)*dy_tr;
        Vsurvey(i_date).VolVE_tr(i_transect)=Vsurvey(i_date).VolVE(i_transect)*dy_tr;
        
    end

    Vsurvey(i_date).VtotalA=sum(Vsurvey(i_date).VolVA_tr);
    Vsurvey(i_date).VtotalB=sum(Vsurvey(i_date).VolVB_tr);
    Vsurvey(i_date).VtotalC=sum(Vsurvey(i_date).VolVC_tr);
    Vsurvey(i_date).VtotalD=sum(Vsurvey(i_date).VolVD_tr);
    Vsurvey(i_date).VtotalE=sum(Vsurvey(i_date).VolVE_tr);

    VtotalA(i_date)=sum(Vsurvey(i_date).VolVA_tr);
    VtotalB(i_date)=sum(Vsurvey(i_date).VolVB_tr);
    VtotalC(i_date)=sum(Vsurvey(i_date).VolVC_tr);
    VtotalD(i_date)=sum(Vsurvey(i_date).VolVD_tr);
    VtotalE(i_date)=sum(Vsurvey(i_date).VolVE_tr);
end












%% Figure Volume changes wrt to first survey in boxes
% Planfom boxes!
figure

plot(datenum(dates(pairs(1,1),:)),0,'or') % zero point
hold on;
bar(datenum(dates(pairs(:,1),:)),(Vtotal-Vtotal(1))/1745)
plot(datenum(dates(pairs(:,1),:)),(VtotalA-VtotalA(1))/1745,'*--g')
plot(datenum(dates(pairs(:,1),:)),(VtotalB-VtotalB(1))/1745,'o--m')
plot(datenum(dates(pairs(:,1),:)),(VtotalC-VtotalC(1))/1745,'+--k')
plot(datenum(dates(pairs(:,1),:)),(VtotalD-VtotalD(1))/1745,'v--b')
plot(datenum(dates(pairs(:,1),:)),(VtotalE-VtotalE(1))/1745,'^--r')
ylim([-250 200])
xlim([ datenum(dates(pairs(1,1),:)) datenum(dates(pairs(end,1),:)) ])
datetick('x')
grid on
ylabel('\Delta V per m alongshore wrt June 2009')
legend('First survey','Total area','Area d A','Area d B','Area d C','Area d D','Area d E','location','northwest')
set(gcf,'PaperPosition',[0  0  25 20])
hold off

title('Cumulative changes since June 2009. Horizontal slices')

if savefigs==1
    filename=[ 'Volume control area y=  ' num2str(B_south) '_' num2str(B_north) ', x=  ' num2str(B_ons_all) '_' num2str(B_offs_all)]
  %  print ([ Mpath '\VolumeAnalysis\images\' filename '.jpg'],'-djpeg')
end

%% Figure first and last topo and areas of interest
figure
% alongshore boundaries
B_north  = 1745;
B_south  = 0;

subplot(3,2,[1 3])
i_date=1;
pcolor(Vsurvey(i_date).transects.y,Vsurvey(i_date).transects.cr_sh,Vsurvey(i_date).transects.elevation);  % plot first survey
shading interp
hold on;
[C,h] = contour(Vsurvey(i_date).transects.y,Vsurvey(i_date).transects.cr_sh,Vsurvey(i_date).transects.elevation,[-9:1:5],'--k');
clabel(C,h,'LabelSpacing',72,'FontSize',8,'Color','k','Rotation',0) ;

rectangle('Position', [B_south B_ons_all (B_north-B_south) (B_offs_all-B_ons_all)],'LineWidth',2,'LineStyle','--')
rectangle('Position', [B_south B_ons_1+2 (B_north-B_south) (B_offs_1-B_ons_1)-2],'EdgeColor','g')
rectangle('Position', [B_south B_ons_2 (B_north-B_south) (B_offs_2-B_ons_2)-5],'EdgeColor','m')
rectangle('Position', [B_south B_ons_3 (B_north-B_south) (B_offs_3-B_ons_3)-5],'EdgeColor','k')
rectangle('Position', [B_south B_ons_4 (B_north-B_south) (B_offs_4-B_ons_4)-5],'EdgeColor','b')
rectangle('Position', [B_south B_ons_5 (B_north-B_south) (B_offs_5-B_ons_5)-5],'EdgeColor','r')
text(B_south+10, B_offs_1-30,'Area I','FontAngle','italic','color','w')
text(B_south+10, B_offs_2-30,'Area II','FontAngle','italic','color','w')
text(B_south+10, B_offs_3-30,'Area III','FontAngle','italic','color','k')
text(B_south+10, B_offs_4-30,'Area IV','FontAngle','italic','color','k')
text(B_south+10, B_offs_5-30,'Area V','FontAngle','italic','color','w')
axis equal
ylim([-10 1210])
caxis([-10 4])
xlabel ('along shore [m]')
ylabel ('cross shore [m]')
title ([dates(pairs(1,1),[1:2 6:end]) ' survey' ])

subplot(3,2,[2 4 ])
i_date=length(Vsurvey);
pcolor(Vsurvey(i_date).transects.y,Vsurvey(i_date).transects.cr_sh,Vsurvey(i_date).transects.elevation);  % plot latest survey
shading interp
hold on;
[C,h] = contour(Vsurvey(i_date).transects.y,Vsurvey(i_date).transects.cr_sh,Vsurvey(i_date).transects.elevation,[-9:1:5],'--k');
clabel(C,h,'LabelSpacing',72,'FontSize',8,'Color','k','Rotation',0) ;
rectangle('Position', [B_south B_ons_all (B_north-B_south) (B_offs_all-B_ons_all)],'LineWidth',2,'LineStyle','--')
rectangle('Position', [B_south B_ons_1+2 (B_north-B_south) (B_offs_1-B_ons_1)-2],'EdgeColor','y')
rectangle('Position', [B_south B_ons_2 (B_north-B_south) (B_offs_2-B_ons_2)-5],'EdgeColor','r')
rectangle('Position', [B_south B_ons_3 (B_north-B_south) (B_offs_3-B_ons_3)-5],'EdgeColor','m')
rectangle('Position', [B_south B_ons_4 (B_north-B_south) (B_offs_4-B_ons_4)-5],'EdgeColor','b')
rectangle('Position', [B_south B_ons_5 (B_north-B_south) (B_offs_5-B_ons_5)-5],'EdgeColor','g')
text(B_south+10, B_offs_1-30,'Area I','FontAngle','italic','color','w')
text(B_south+10, B_offs_2-30,'Area II','FontAngle','italic','color','w')
text(B_south+10, B_offs_3-30,'Area III','FontAngle','italic','color','k')
text(B_south+10, B_offs_4-30,'Area IV','FontAngle','italic','color','k')
text(B_south+10, B_offs_5-30,'Area V','FontAngle','italic','color','w')
axis equal
ylim([-50 1300])
caxis([-10 4])
xlabel ('along shore [m]')
title ([dates(pairs(end,1),[1:2 6:end]) ' survey' ])

subplot(3,2,[5 6])
caxis([-10 4])
colorbar('location', 'north')
axis off

set(gcf,'PaperPosition',[0  0  25 20])
print Overview_firstandlastsurveyPlanformboxes.png -dpng



%% relate to wave power
load('waveheightanalysis/Pperiod.mat')

%Pperiod(:,4)=Pperiod(:,1)./(Pperiod(:,3)-Pperiod(:,2));

% timeserie
figure(34)
hold on;
for ii=2:length(m_p)
    bar(datenum(dates(pairs(ii,1),:)),(m_p(ii)-m_p(ii-1))*(B_offs_all-B_ons_all) )
    plot(datenum(dates(pairs(ii,1),:)),(m1_p(ii)-m1_p(ii-1))*(B_offs_1-B_ons_1),'*y')
    plot(datenum(dates(pairs(ii,1),:)),(m2_p(ii)-m2_p(ii-1))*(B_offs_2-B_ons_2),'om')
    plot(datenum(dates(pairs(ii,1),:)),(m3_p(ii)-m3_p(ii-1))*(B_offs_3-B_ons_3),'vr')
    plot(datenum(dates(pairs(ii,1),:)),(m4_p(ii)-m4_p(ii-1))*(B_offs_4-B_ons_4),'^b')
    plot(datenum(dates(pairs(ii,1),:)),(m5_p(ii)-m5_p(ii-1))*(B_offs_5-B_ons_5),'xg')
end

plot( Pperiod(2:end,3), Pperiod(2:end,1)/10^10)
plot( Pperiod(2:end,3), Pperiod(2:end,4)/10^9)

ylim([-120 120])
xlim([ datenum(dates(pairs(1,1),:)) datenum(dates(pairs(end,1),:)) ])
datetick('x')
grid on
ylabel('\Delta V per m alongshore wrt June 2009')
legend('Total area','Area I','Area II','Area III','Area IV','Area V','location','northwest')

set(gcf,'PaperPosition',[0  0  25 20])

% Scatter
figure(35)
hold on;
for ii=2:length(m_p)
    scatter(Pperiod(ii,4),(m_p(ii)-m_p(ii-1))*(B_offs_all-B_ons_all),'.k' )
    scatter(Pperiod(ii,4),(m1_p(ii)-m1_p(ii-1))*(B_offs_1-B_ons_1),'*y')
    scatter(Pperiod(ii,4),(m2_p(ii)-m2_p(ii-1))*(B_offs_2-B_ons_2),'om')
    scatter(Pperiod(ii,4),(m3_p(ii)-m3_p(ii-1))*(B_offs_3-B_ons_3),'vr')
    scatter(Pperiod(ii,4),(m4_p(ii)-m4_p(ii-1))*(B_offs_4-B_ons_4),'^b')
    scatter(Pperiod(ii,4),(m5_p(ii)-m5_p(ii-1))*(B_offs_5-B_ons_5),'xg')
end

%  plot( Pperiod(2:end,3), Pperiod(2:end,1)/10^10)


%  ylim([-120 120])
%  xlim([ datenum(dates(pairs(1,1),:)) datenum(dates(pairs(end,1),:)) ])
%  datetick('x')
grid on
ylabel('\Delta V per m alongshore wrt June 2009')
legend('Total area','Area I','Area II','Area III','Area IV','Area V','location','northwest')

set(gcf,'PaperPosition',[0  0  25 20])





%% Figure Volume changes in boxes
figure(32)

for pair_nr=2:19;
    dV_all(pair_nr)=(m_p(pair_nr)-m_p(pair_nr-1))*(B_offs_all-B_ons_all);
    dV_1(pair_nr)=(m1_p(pair_nr)-m1_p(pair_nr-1))*(B_offs_1-B_ons_1);
    dV_2(pair_nr)=(m2_p(pair_nr)-m2_p(pair_nr-1))*(B_offs_2-B_ons_2);
    dV_3(pair_nr)=(m3_p(pair_nr)-m3_p(pair_nr-1))*(B_offs_3-B_ons_3);
    dV_4(pair_nr)=(m4_p(pair_nr)-m4_p(pair_nr-1))*(B_offs_4-B_ons_4);
    dV_5(pair_nr)=(m5_p(pair_nr)-m5_p(pair_nr-1))*(B_offs_5-B_ons_5);
end

bar(datenum(dates(pairs(:,1),:)),dV_all)
hold on;
plot(datenum(dates(pairs(:,1),:)),dV_1,'*--y')
plot(datenum(dates(pairs(:,1),:)),dV_2,'m')
plot(datenum(dates(pairs(:,1),:)),dV_3,'r')
plot(datenum(dates(pairs(:,1),:)),dV_4,'b')
plot(datenum(dates(pairs(:,1),:)),dV_5,'g')
ylim([-250 200])
xlim([ datenum(dates(pairs(1,1),:)) datenum(dates(pairs(end,1),:)) ])
datetick('x')

grid on
ylabel('\Delta V per m alongshore wrt previous survey')
legend('Total area','Area I','Area II','Area III','Area IV','Area V','location','northwest')

set(gcf,'PaperPosition',[0  0  25 20])
print Vol_change_total.png -dpng

%% Figure mean elevation changes in boxes
figure

bar(datenum(dates(pairs(:,1),:)),m_p-m_p(1))
hold on;
plot(datenum(dates(pairs(:,1),:)),m1_p-m1_p(1),'y')
plot(datenum(dates(pairs(:,1),:)),m2_p-m2_p(1),'m')
plot(datenum(dates(pairs(:,1),:)),m3_p-m3_p(1),'r')
plot(datenum(dates(pairs(:,1),:)),m4_p-m4_p(1),'b')
plot(datenum(dates(pairs(:,1),:)),m5_p-m5_p(1),'g')
ylim([-1.2 1])
xlim([ datenum(dates(pairs(1,1),:)) datenum(dates(pairs(end,1),:)) ])
datetick('x')
grid on
ylabel('mean elevation change')
legend('Total area','Area I','Area II','Area III','Area IV','Area V','location','northwest')

set(gcf,'PaperPosition',[0  0  25 20])
print Elevation_Areas.png -dpng


%% figures all surveys

%% Figure first and last topo and areas of interest
if 0
    figure(322)
    for ipair=1:length(pairs)
        pcolor(xi,yi,altitude_i(:,:,ipair));  % plot first survey
        shading interp
        hold on;
        [C,h] = contour(xi,yi,altitude_i(:,:,ipair),[-9:1:5],'--k');
        clabel(C,h,'LabelSpacing',72,'FontSize',8,'Color','k','Rotation',0) ;
        
        rectangle('Position', [B_south B_ons_all (B_north-B_south) (B_offs_all-B_ons_all)],'LineWidth',2,'LineStyle','--')
        rectangle('Position', [B_south B_ons_1+2 (B_north-B_south) (B_offs_1-B_ons_1)-2],'EdgeColor','y')
        rectangle('Position', [B_south B_ons_2 (B_north-B_south) (B_offs_2-B_ons_2)-5],'EdgeColor','r')
        rectangle('Position', [B_south B_ons_3 (B_north-B_south) (B_offs_3-B_ons_3)-5],'EdgeColor','m')
        rectangle('Position', [B_south B_ons_4 (B_north-B_south) (B_offs_4-B_ons_4)-5],'EdgeColor','b')
        rectangle('Position', [B_south B_ons_5 (B_north-B_south) (B_offs_5-B_ons_5)-5],'EdgeColor','g')
        text(B_south+10, B_offs_1-30,'Area I','FontAngle','italic','color','w')
        text(B_south+10, B_offs_2-30,'Area II','FontAngle','italic','color','w')
        text(B_south+10, B_offs_3-30,'Area III','FontAngle','italic','color','k')
        text(B_south+10, B_offs_4-30,'Area IV','FontAngle','italic','color','k')
        text(B_south+10, B_offs_5-30,'Area V','FontAngle','italic','color','w')
        axis equal
        ylim([-50 1300])
        caxis([-10 4])
        xlabel ('along shore [m]')
        ylabel ('cross shore [m]')
        title ([dates(pairs(ipair,1),[4:6 1:3 7:end]) ' Survey' ])
        colorbar
        
        set(gcf,'Renderer','zbuffer')
        set(gcf,'PaperPosition',[0  0  25 20])
        filenameComb=['CombinedSurvey' dates(pairs(ipair,1),[7:end 1:2 4:5 ]) '.png']
        print (filenameComb, '-dpng')
        hold off
    end
end


