% 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_000000_rst.nc'];
ncmap2   = [datadir caseid '_20151101_001000_rst.nc'];
ncmap3   = [datadir caseid '_20151101_002000_rst.nc'];
ncmap   = [datadir caseid '_map.nc'];

info=ncinfo(ncmap);
info1=ncinfo(ncmap1);
info2=ncinfo(ncmap2);
info3=ncinfo(ncmap3);

%ncdisp(ncmap);
time=ncread(ncmap,'time');
%1
FlowElem_xbnd1    = ncread(ncmap1,'FlowElem_xbnd');
FlowElem_ybnd1    = ncread(ncmap1,'FlowElem_ybnd');
bl_bnd1 = ncread(ncmap1,'bl_bnd'); 
bl_bnd1 = bl_bnd1(:,end);
id_1D2D1 = find(FlowElem_ybnd1 <-4); % find which results below to the 1D2D grid (here y<-4 is used)
bl_bnd1 =  bl_bnd1(id_1D2D1,:);
FlowElem_xbnd1 = FlowElem_xbnd1(id_1D2D1,:);
FlowElem_ybnd1 = FlowElem_ybnd1(id_1D2D1,:);
%2
FlowElem_xbnd2    = ncread(ncmap2,'FlowElem_xbnd');
FlowElem_ybnd2    = ncread(ncmap2,'FlowElem_ybnd');
bl_bnd2 = ncread(ncmap2,'bl_bnd'); 
bl_bnd2 = bl_bnd2(:,end);
id_1D2D2 = find(FlowElem_ybnd2 <-4); % find which results below to the 1D2D grid (here y<-4 is used)
bl_bnd2 =  bl_bnd2(id_1D2D2,:);
FlowElem_xbnd2 = FlowElem_xbnd2(id_1D2D2,:);
FlowElem_ybnd2 = FlowElem_ybnd2(id_1D2D2,:);
%3
FlowElem_xbnd3    = ncread(ncmap3,'FlowElem_xbnd');
FlowElem_ybnd3    = ncread(ncmap3,'FlowElem_ybnd');
bl_bnd3 = ncread(ncmap3,'bl_bnd'); 
bl_bnd3 = bl_bnd3(:,end);

id_1D2D3 = find(FlowElem_ybnd3 <-4); % find which results below to the 1D2D grid (here y<-4 is used)
bl_bnd3 =  bl_bnd3(id_1D2D3,:);
FlowElem_xbnd3 = FlowElem_xbnd3(id_1D2D3,:);
FlowElem_ybnd3 = FlowElem_ybnd3(id_1D2D3,:);

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







% 
% 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);
plot(FlowElem_xbnd1,bl_bnd1,'--gs',FlowElem_xbnd2,bl_bnd2,'b-o',FlowElem_xbnd3,bl_bnd3,'r*');
grid on
%pp=plot(xbel,hbel-d1,'-dk'); %set(pp,'LineWidth',2);
legend('t0','t10','t20');
% 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])
xlabel(' x-distance along the us-boundary [m]')
%xlabel(' y-distance')
ylabel(' bed level at the boundary [m]')

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