% Program to plot velocity magnitudes
% Andries Paarlberg
% paarlberg@hkv.nl 
% 24-11-2017

% vard = vardata;

varname  = 'vel. magnitude';
varunit  = '(m s^{-1})';
axlims   = [0 30 0.7 0.9];
%axlims = [0 30 -Inf +Inf]
legpos   = 'best';

% 1D-data (no manipulation needed)
xcor_1D  = mesh1d_node_x;
vard_1D  = 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_ucmag(:,id_time));
vard_2D   = squeeze(mesh2d_ucmag(:,id_time));

%% 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( mesh2d_face_y>1.30 & mesh2d_face_y<1.40 );
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;

%% plotting

%% make figure
figure
subplot(2,1,1)
hold on; grid on; box on;
plot( xcor_1D , vard_1D , '.')
plot( xcor_F2D, vard_F2D, 'x')
title([{[strrep(tit_text,'_','\_') ' ' varname]},' comparison 1D and 2D'],'fontsize',10)
xlabel('x (m)','fontsize',10)
ylabel([{varname} ;{varunit}],'fontsize',10)
legend('1D','2D','location',legpos)
axis(axlims)

% subplot(3,1,2)
% hold on; grid on; box on;
% %plot( xcor_Q2D, vard_Q2D, '.')
% plot( xcor_F2D, vard_F2D, 'x')
% title([{[tit_text varname]},' 2D comparison '],'fontsize',10)
% xlabel('x (m)','fontsize',10)
% ylabel([{varname} ;{varunit}],'fontsize',10)
% legend('2D','quasi2D','location',legpos)
% axis(axlims)

subplot(2,1,2)
hold on; grid on; box on;
%plot( xcor_Q2D, vard_F2D - vard_Q2D, '.')
plot( xcor_1D , vard_1D  - vard_F2D, '.')
title([{[strrep(tit_text,'_','\_') ' ' varname]},' differences '],'fontsize',10)
xlabel('x (m)','fontsize',10)
ylabel([{['diff. ' varname]} ;{varunit}],'fontsize',10)
xlim(axlims(1:2))
legend('1D minus 2D','location',legpos)
%axis(axlims)

% overal title 

% store figure
x0=11;
y0=0;
if f_sav ==1;
    width=16.5;
    height=20.5;
    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 = [caseid '_' varname ,'_map', num2str(id_time)];
    figname = strrep(figname,' ','');
    figname = strrep(figname,'.','');
    print([figdir,'/pdf/',figname,'.pdf'],'-dpdf','-r300')
    print([figdir,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
