% Compare water depth for steady uniform flow with zero bed slope
% Andries Paarlberg/ modified by A. Omer
% 21-11-2017

clear all;
close all;

if isempty('vs_ini')
    run('d:\openearth\matlab\oetsettings.m')
end

caseid  = 'str01';
datadir = '../dflowfmoutput/';
ncmap   = [datadir caseid '_map.nc'];

info=ncinfo(ncmap);
%ncdisp(ncmap);
time=ncread(ncmap,'time');

mesh2d_face_x     = ncread(ncmap,'mesh2d_face_x');
mesh2d_face_y     = ncread(ncmap,'mesh2d_face_y');

mesh2d_ucmag = ncread(ncmap,'mesh2d_ucmag'); 
mesh2d_ucmag = mesh2d_ucmag(:,end);

mesh1d_node_x     = ncread(ncmap,'mesh1d_node_x');
mesh1d_node_y     = ncread(ncmap,'mesh1d_node_y');

mesh1d_ucmag = ncread(ncmap,'mesh1d_ucmag');
mesh1d_ucmag = mesh1d_ucmag(:,end);

[mesh1d_node_x,I] = sortrows(mesh1d_node_x);
mesh1d_ucmag = mesh1d_ucmag(I);
clear I;

d1 = 0.5;   % water depth at downstream end

y0 = -5;
y1 = 11.;
dy = 0.16;
y  = [y0+dy/2:dy:y1-dy/2]';

y(1)   = y(1)  +1e-8;
y(end) = y(end)-1e-8;

xc = -4.2;                     % y-coordinate centerline
d_2Dfull = griddata(mesh2d_face_y,mesh2d_face_x,mesh2d_ucmag,y,repmat(xc,length(y),1));

d_1D     = mesh1d_ucmag;



% %% belanger
% ybel0 = -5;
% L  = 16;
% ib = 0;
% C  = 45;
% B  = 0.08;
% hb = 0.67;   % waterdiepte op benedenstroomse rand
% 
% %Q = B * hb * C * sqrt( hb * ib );
% Q = 0.08;     % input instead of steady flow equivalent
% N = 100;      % for Belanger-computation
% 
% %function [x,h] = belanger(Q,C,L,B,ib,hb,N);
% [ybel,hbel] = belanger(Q,C,L,B,ib,hb,N);
% ybel=ybel';
% hbel=hbel';
% ybel=ybel+ybel0;

%figure;
%plot(x,h);

% mean(hbel)

figure;
subplot(2,1,1);
hold on; grid on; box on;
pp=plot(y,d_2Dfull,'-xr');
pp=plot(mesh1d_node_y,d_1D,'-dg');
%pp=plot(xbel,hbel-d1,'-dk'); %set(pp,'LineWidth',2);
legend('2D','1D','Location', 'southeast');
ylabel(' Flow element cent.velocity magnitude [m/s]')
xlim([-5 11])

subplot(2,1,2);
hold on; grid on; box on;
pp=plot(y,[flipud(d_1D) - d_2Dfull],'-xr');
%pp=plot(mesh1d_node_x,d_1D-hbel(2:end-1),'-dg');
%legend('2D','1D','location','southeast');
ylabel('velocity difference (1D-2D)[m/s]')
xlim([-5 11])
xlabel('distance along y-direction [m]')
print(gcf,'-dpng',[caseid '.png'])
