% 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/';
ncmap1   = [datadir caseid '_20151101_002000_rst.nc'];
ncmap   = [datadir caseid '_map.nc'];

info=ncinfo(ncmap);
info1=ncinfo(ncmap1);

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

FlowElem_xbnd    = ncread(ncmap1,'FlowElem_xbnd');
FlowElem_ybnd    = ncread(ncmap1,'FlowElem_ybnd');

bl_bnd = ncread(ncmap1,'bl_bnd'); 
bl_bnd = bl_bnd(:,end);


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

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

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

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

[mesh1d_node_x,I] = sortrows(mesh1d_node_x);
mesh1d_mor_bl = mesh1d_mor_bl(I);

d1 = 0.5;   % water depth at downstream end

y0 = -5;
x1 = 11.;
dy = 0.16;
y  = [y0+dy/2:dy:x1-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_mor_bl,y,repmat(xc,length(y),1));

d_1D     = mesh1d_mor_bl;


%% 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)
x1(1:100)= xc;
%  figure;
% subplot(2,1,1);
% hold on;
%grid on
% pp=plot3(x,y1,d_2Dfull-d1,'-xr');
% pp=plot3(mesh1d_node_x,y1, d_1D-d1,'-dg');
% pp=scatter3(FlowElem_xbnd, FlowElem_ybnd,bl_bnd);
plot3(y,x1,d_2Dfull,'--gs',mesh1d_node_y,x1,d_1D,'b-o',FlowElem_ybnd, FlowElem_xbnd,bl_bnd,'c*');
grid on
%pp=plot(xbel,hbel-d1,'-dk'); %set(pp,'LineWidth',2);
legend('2D','1D','bnd-bl');
% ylabel('Bed level [m]')
xlim([-5 11])

% subplot(2,1,2);
% hold on; grid on; box on;
% pp=plot(x,[d_2Dfull-d_1D],'-xr');
% %pp=plot(mesh1d_node_x,d_1D-hbel(2:end-1),'-dg');
% %legend('2D','1D','location','southeast');
% ylabel('Bed level difference (2D-1D)[m]')
% xlim([0 31])
ylabel(' x-direction')
xlabel(' y-distance')
zlabel(' bed level')

print(gcf,'-dpng',[caseid '.png'])
