% Compare water depth/level for steady uniform flow with zero bed slope
% 
% Andries Paarlberg
% 24-11-2017

clear all;
close all;

if ~exist('vs_ini')
    run('d:\checkouts\openearthtools\matlab\oetsettings.m')
end

caseid  = 't01';
datadir = '../dflowfmoutput/';
fig_dir = '../doc/figures/';
ncmap   = [datadir caseid '_map.nc'];

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

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

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

x = mesh1d_node_x;
y = mesh1d_node_y;
z = mesh1d_waterdepth;

id1 = find( y < 2 );
id2 = find( y > 4 & x < 21 );
id3 = find( y > 4 & x > 24 & x < 41 );
id4 = find( y > 4 & x > 44 );

x1 = x(id1); y1 = y(id1); z1 = z(id1);
x2 = x(id2); y2 = y(id2); z2 = z(id2);
x3 = x(id3); y3 = y(id3); z3 = z(id3);
x4 = x(id4); y4 = y(id4); z4 = z(id4);

format long;
L1 = pathdistance(x1,y1); L1(end)
L2 = pathdistance(x2,y2); L2(end)
L3 = pathdistance(x3,y3); L3(end)
L4 = pathdistance(x4,y4); L4(end)

d1 = 0.35;   % water depth at downstream end

%% belanger
xbel0 = 0;
L  = 30;
ib = 0;
C  = 45;
B  = 0.3;
hb = d1;   % 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);
[xbel,hbel] = belanger(Q,C,L,B,ib,hb,N);
xbel=xbel';
hbel=hbel';
xbel=xbel+xbel0;

%figure;
%plot(x,h);

mean(hbel)

dmax = max([z1;z2;z3;z4;hbel]);

figure;
hold on; grid on; box on;
pp=plot(xbel,hbel,'-k'); set(pp,'LineWidth',2);
plot(L1,z1,'-r')
plot(L2,z2,'-g')
plot(L3,z3,'-b')
plot(L4,z4,'-m')
legend('analytic dx0.50 straight','dx0.50m straight','dx0.25m 90gr bend','dx0.50m 90gr bend','dx1.00m 90gr bend','location','best');
xlabel('x [m]')
ylabel('water depth [m]')
ylim([d1-0.002 dmax+0.002])
xlim([0 31])

%% export fig
x0=11;
y0=0;
width=15;
height=10;
tit_text = '';
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 = ['compare_waterdepth'];
figname = strrep(figname,' ','');
figname = strrep(figname,'.','');
%print([fig_dir,'pdf/',figname,'.pdf'],'-dpdf','-r300')
print([fig_dir,figname,'.png'],'-dpng','-r300')
