% Program to plot velocity magnitudes
% Andries Paarlberg
% paarlberg@hkv.nl 
% 24-11-2017
% adjusted by Stef Boersen
% stef.boersen@rhdhv.com
% 29-11-2017 

% vard = vardata;

varname  = 'debiet';
varunit  = '(m^{3} s^{-1})';
axlims   = [0 30 0.029999999995 0.030000000005];
%axlims = [0 30 -Inf +Inf]
legpos   = 'best';

% 1D-data (no manipulation needed)
xcor_1D  = mesh1d_node_x;
vard_1D  = mesh1d_waterdepth.*mesh1d_ucmag;

% 2D grid
xcor_2D  = mesh2d_face_x;
ycor_2D  = mesh2d_face_y;

% find time id, if last mapfile requested
if id_time == -1
    % find last time on mapfile
    id_time = size(vard_1D,2);
end

% reduce data to requested map on ncfiles
vard_1D   = squeeze(mesh1d_waterdepth(:,id_time).*mesh1d_ucmag(:,id_time)*0.1);
vard_2D   = squeeze(mesh2d_waterdepth(:,id_time).*mesh2d_ucmag(:,id_time)*0.1);

%% adjust the data 2D-data
% do some manipulation on 2D-data to differentiate between full and qausi
% 2D-grid (specific for this case!)

% find elements of Quasi2D grid (2D, 1 cell wide)
% results below to the 1D2D grid (here y=1 is used)
id_Q2D    = find(ycor_2D>1);
xcor_Q2D  = xcor_2D(id_Q2D);
vard_Q2D  = vard_2D(id_Q2D);

% find elements of Full2D grid (2D, 5 cells wide)
% take centerline
id_F2D    = find( ycor_2D>0.21 & ycor_2D<0.29 );
xcor_F2D  = xcor_2D(id_F2D);
vard_F2D  = vard_2D(id_F2D);

% sort based on xcor
% 1D
[~,sort_idx] = sortrows(xcor_1D);
xcor_1D  = xcor_1D(sort_idx);
vard_1D  = vard_1D(sort_idx);
% Quasi2D
[~,sort_idx] = sortrows(xcor_Q2D);
xcor_Q2D  = xcor_Q2D(sort_idx);
vard_Q2D  = vard_Q2D(sort_idx);
% Full2D
[~,sort_idx] = sortrows(xcor_F2D);
xcor_F2D  = xcor_F2D(sort_idx);
vard_F2D  = vard_F2D(sort_idx);
clear sort_idx;

%% Add values from his files
% find timeindex in his file
id_time_cross=find(time_cross==time(id_time));

% select correct indexen
id_1D_cross  = find(y_cor_cross(2,:)>2.4 & y_cor_cross(2,:)< 3.4 );  
id_Q2D_cross = find(y_cor_cross(2,:)>1.0 & y_cor_cross(2,:)< 2.4 );
id_F2D_cross = find(y_cor_cross(2,:)>-0.01 & y_cor_cross(2,:)<0.54 );

% find elements of Full2D grid (2D, 5 cells wide)
% take centerline

xcor_1D_cross  = x_cor_cross(1,id_1D_cross);
xcor_Q2D_cross = x_cor_cross(1,id_Q2D_cross);
xcor_F2D_cross = x_cor_cross(1,id_F2D_cross);

var_1D_cross  =dis_cross(id_1D_cross,id_time_cross);
var_Q2D_cross =dis_cross(id_Q2D_cross,id_time_cross);
var_F2D_cross =dis_cross(id_F2D_cross,id_time_cross)/5; %deviation by 5 required because of number of gridcells
%% plotting

%% make figure
figure
subplot(3,1,1)
hold on; grid on; box on;
plot( xcor_1D , vard_1D , 'r.')
plot( xcor_Q2D, vard_Q2D, 'b.')
plot( xcor_Q2D_cross,var_Q2D_cross,'bO')
plot( xcor_1D_cross,var_1D_cross,'rx')
title([{[tit_text varname]},' 1D comparison '],'fontsize',10)
xlabel('x (m)','fontsize',10)
ylabel([{varname} ;{varunit}],'fontsize',10)
legend('1D map','quasi2D map','1D his','quasi2D his','location',legpos)
axis(axlims)

subplot(3,1,2)
hold on; grid on; box on;
plot( xcor_Q2D, vard_Q2D, 'b.')
plot( xcor_F2D, vard_F2D, 'r.')
plot( xcor_Q2D_cross,var_Q2D_cross,'bO')
plot( xcor_F2D_cross,var_F2D_cross,'rx')
title([{[tit_text varname]},' 2D comparison '],'fontsize',10)
xlabel('x (m)','fontsize',10)
ylabel([{varname} ;{varunit}],'fontsize',10)
legend('2D map','quasi2D map','2D his','quasi2D his','location',legpos)
axis(axlims)

subplot(3,1,3)
hold on; grid on; box on;
plot( xcor_Q2D, vard_F2D - vard_Q2D, 'b.')
plot( xcor_1D , vard_1D  - vard_Q2D, 'r.')
plot( xcor_Q2D_cross , var_F2D_cross - var_Q2D_cross, 'bO')
plot( xcor_Q2D_cross , var_1D_cross  - var_Q2D_cross, 'rx')
title([{[tit_text varname]},' differences '],'fontsize',10)
xlabel('x (m)','fontsize',10)
ylabel([{['diff. ' varname]} ;{varunit}],'fontsize',10)
xlim(axlims(1:2))
legend('2D minus Quasi2D map','1D minus Quasi2D map','2D minus Quasi2D his','1D minus Quasi2D his','location',legpos)
ylim([-5e-16 5e-16])
%axis(axlims)

% overal title 

% store figure
x0=11;
y0=0;
if f_sav ==1;
    width=16.5;
    height=24.7;
    tit_2=tit_text(find(~isspace(tit_text))); %remove whitespaces from name 
    set(gcf,'units','centimeters','position',[x0,y0,width,height]);
    set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[width, height]);  
    figname = [varname ,'_map', num2str(id_time)];
    figname = strrep(figname,' ','');
    figname = strrep(figname,'.','');
    print([fig_dir,figname,'.pdf'],'-dpdf','-r300')
    print([fig_dir,figname,'.png'],'-dpng','-r300')
end

%% opruimen
if f_debug ~=1
clearvars id_1D2D bedload_1D2D x_cor1D2D y_cor1D2D id_2D bedload_2D_full ...
            x_cor2D_full y_cor2D_full
end
