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

clear all;
close all;

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

caseid  = 'belanger_AP_rst';
datadir = ['../DFM_OUTPUT_' caseid '/'];
%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_waterdepth = ncread(ncmap,'mesh2d_waterdepth'); 
mesh2d_waterdepth = mesh2d_waterdepth(:,end);

mesh2d_waterlevel = ncread(ncmap,'mesh2d_s1'); 
mesh2d_waterlevel = mesh2d_waterlevel(:,end);

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

%mesh1d_waterdepth = ncread(ncmap,'mesh1d_waterdepth');
%size(mesh1d_waterdepth,2)
%mesh1d_waterdepth = mesh1d_waterdepth(:,end);

%d_1D     = mesh1d_waterdepth;

%[mesh1d_node_x,I] = sortrows(mesh1d_node_x);
%mesh1d_waterdepth = mesh1d_waterdepth(I);
%clear I;

[mesh2d_face_x,I] = sortrows(mesh2d_face_x);
%mesh2d_waterdepth = mesh2d_waterdepth(I);
mesh2d_waterlevel = mesh2d_waterlevel(I);
clear I;

BEDLEV = interp1([0 100000],[0 -10],mesh2d_face_x);

d_2D     = mesh2d_waterlevel - BEDLEV;

%hb = 10 -0.12600 ;   % water depth at downstream end
hb = -0.12600 - (-10-0.025) ;   % water depth at downstream end
xb = 100250;  % locatie waar deze diepte geldt (halve roostercel verschoven)

%% belanger
ib = 1e-4;
C  = 60;
B  = 20;
L  = 100500; %.+500;

%Q = B * hb * C * sqrt( hb * ib );
Q = 600;      % input instead of steady flow equivalent
N = 201;      % for Belanger-computation

%function [x,h] = belanger(Q,C,L,B,ib,hb,N);
[xbel,hbel] = belanger_dx1m(Q,C,L,B,ib,xb,hb,N);
xbel=xbel';
hbel=hbel';

%xbel=xbel+xbel0;

%figure;
%plot(x,h);

%mean(hbel)

figure;
subplot(2,1,1);
hold on; grid on; box on;
pp=plot(mesh2d_face_x,d_2D,'-xr');
%pp=plot(mesh1d_node_x,d_1D,'-dg');
pp=plot(xbel,hbel,'-dk'); %set(pp,'LineWidth',2);
legend('2D','analytic','location','best');
%legend('1D','analytic');
ylabel('water depth [m]')
xlim([ -500 1e5+500 ])

subplot(2,1,2);
hold on; grid on; box on;
pp=plot(mesh2d_face_x,d_2D-hbel(2:end-1),'-xr');
%pp=plot(mesh1d_node_x,d_1D-hbel(2:end-1),'-dg');
legend('2D','location','best');
ylabel('water depth diff wrt analytic[m]')
xlim([ -500 1e5+500 ])
%xlim([0 31])

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