%% Script to investigate the error in the bathymetry due to a poor sounding
%% velocity, and correct it.

% approach: 
% 1. For each survey, calculate the mean volume in the offshore area where no bed level changes are to be expected (depth > 7 m)
% 2. correct it by suggesting a multiplication factor.


%%  boundaries Area to be investigated

B_top    = -7;  % top boundary 
B_bottom = -8.5;   % bottom boundary



%% locate vlugtenburg directory on this computer
curdir=pwd;
ind_v=findstr(curdir, 'vlugtenburg');
pathvlugtenburg=(curdir(1:ind_v+11));  % gives back cometing like F:\check_outs\tudelft\vlugtenburg\

%%  load survey data 
close all
NetCDFfilename=[pathvlugtenburg 'Vlugtenburg.nc'];
%nc_dump(NetCDFfilename);  % give a quick overview of the contents

dates=nc_varget(NetCDFfilename,'time');     % index -1 because netcdf starts counting at zero.. Weird Yanks.
alt=nc_varget(NetCDFfilename,'altitude');  % Format is surveydate,linenr,xloc
surveytype=nc_varget(NetCDFfilename,'survey_type');
x=nc_varget(NetCDFfilename,'cross_shore');
survey_path=nc_varget(NetCDFfilename,'survey_path');

Surveys=list_surveys;
indjet=Surveys.indjet; indwalk=Surveys.indwalk; pairs=Surveys.Survey_pairs;  % rewrite to local variables... bit easier later on


%% Show a plot of the last Survey with the box


figure
pcolor(x,squeeze(alt(indjet(end),:,:)));
hold all
%legend((dates(indjet,:)))
% plot(x,x*0+B_top,'--k','linewidth',2);
% plot(x,x*0+B_bottom,'--k','linewidth',2);



%% Calculate Volumes 

alt_bound=nan(size(squeeze(alt(:,iline,:))));

for i_d=1:length(indjet) % for all surveys
    % for each survey do this:
    
xi_bnd=find(squeeze(alt(indjet(i_d),iline,:)) >= B_bottom & squeeze(alt(indjet(i_d),iline,:)) <= B_top) % find the x locations of points where the z is between the boundaries
alt_bound(indjet(i_d),xi_bnd)=squeeze(alt(indjet(i_d),iline,xi_bnd))-B_bottom; % store the profile between the boundaries
dx=x(2)-x(1); %% gridsize in x
V(indjet(i_d))=nansum(alt_bound(indjet(i_d),xi_bnd)*dx)+x(xi_bnd(1))*(abs(B_top-B_bottom)) ; % volume per meter alongshore. (-1 because z values are negative)
end


%% Plot the Transects

figure
plot(x,squeeze(alt(indjet,iline,:)));
hold all
legend((dates(indjet,:)))
plot(x,x*0+B_top,'--k','linewidth',2);
plot(x,x*0+B_bottom,'--k','linewidth',2);


%%
tmin=datenum('01.05.2009', 'dd.mm.yyyy');
tmax=datenum('01.05.2010', 'dd.mm.yyyy');

figure
% subplot top: wave conditions
subplot(6,1,[1:2])
%plot(Hs(:,1),Hs(:,2),'r');
hold on
bar(datenum(dates(indjet,:)),V(indjet)-V(indjet(1)))
legend('Volume change wrt first survey')
for idate=1:length(dates)
    if surveytype(idate,1) == 'J'
        plot([ datenum(dates(idate,:)) datenum(dates(idate,:))], [0 max(V(indjet)-V(indjet(1)))],'-.k')
    end
end
xlim([tmin tmax])
ylabel('m^3 m^{-1}')
datetick('x',20, 'keeplimits')

title ('Onshore migration nourishment Vlugtenburg')


% subplot bottom : timestack surveyline
opt2=1;  %% choose between a timestack of the depths (opt2=1) or a timestack of the sed erosion (opt2=2)

subplot(6,1,[3:6])

if opt2==1
pcolor(datenum(dates(indjet,:)),x,squeeze(alt(indjet,iline,:))'); 
hold on;
pcolor(datenum(dates(indwalk,:)),x,squeeze(alt(indwalk,iline,:))');
shading flat
caxis([-8 1]); %colorbar
% Show line top nourishment
%plot(datenum(dates(indjet,:)),x(ind_nour),'w','linewidth',2)
end

if opt2==2
SedEr_j=diff(squeeze(alt(indjet,iline,:)));      
SedEr_w=diff(squeeze(alt(indwalk,iline,:)));    
pcolor(datenum(dates(indjet(2:end),:)),x,SedEr_j'); 
hold on;
pcolor(datenum(dates(indwalk(2:end),:)),x,SedEr_w');
shading flat
caxis([-1 1]);
cmap=colormap;
cmap(28:35,:)=1;
colormap(cmap)
%colorbar
end


% plot lines in the timestack where the surveys happened
for idate=1:length(dates)
    if surveytype(idate,1) == 'J'
        plot([ datenum(dates(idate,:)) datenum(dates(idate,:))], [0 800],'-.k')
    end
end

ylim([100 800])
xlim([tmin tmax])
datetick('x',20, 'keeplimits')
colorbar('location','West')
xlabel('date')
ylabel('crossshore dist')
title(['Timestack profile evolution line ' num2str(iline)])
finame=(['images\Onshoremigrationnourishment' ]);
print (finame,'-dmeta')

