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

clc;
clear all;
close all;

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

figdir  = '../../c10_straight_channel_1Dvs2D_with_bedslope_t03_Wup_blev1/doc/figures/';

% caseid  = 't02_WLup_bedlevtype1';
% caseid2  = 't02\_WLup\_bedlevtype1';

% caseid  = 't02_WLup_bedlevtype3';
% caseid2  = 't02\_WLup\_bedlevtype3';

% caseid  = 't02_Qup_bedlevtype1';
% caseid2  = 't02\_Qup\_bedlevtype1';

caseid  = 't02_Qup_bedlevtype3';
caseid2  = 't02\_Qup\_bedlevtype3';

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

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

%% read data
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-0);

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-0);
mesh1d_ucmag = ncread(ncmap,'mesh1d_ucmag'); 
%velocity magnitude 1D-grid
%location = 'node'
%SIZE:  99x4
mesh2d_ucmag = ncread(ncmap,'mesh2d_ucmag'); 
%velocity magnitude 2D-grid
%location = 'face' (ofwel de 2D-tegel, daarom:  'Flow element center velocity magnitude')
%SIZE: 594x4

%% plot velocity

f_debug = 1;
f_sav = 1;
tit_text= caseid2;
id_time = -1;
P03_plotoption3;

%% manipulate water level data

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

%d1 = 0.402;
d1 = 0;

x0 = 0.3;
x1 = 30.;
dx = 0.3;
x  = [x0+dx/2:dx:x1-dx/2]';
x(1)   = x(1)  +1e-8;
x(end) = x(end)-1e-8;

yc = 0.25;                     % y-coordinate centerline
d_2Dfull = griddata(mesh2d_face_x,mesh2d_face_y,mesh2d_waterdepth,x,repmat(yc,length(x),1));
yc = 1.75;                     % y-coordinate centerline
d_2D1cel = griddata(mesh2d_face_x,mesh2d_face_y,mesh2d_waterdepth,x,repmat(yc,length(x),1));

d_1D     = mesh1d_waterdepth;

%% belanger
%hb = 0.35 ;    % water depth at downstream end
xb = 30+0.15;  % locatie waar deze diepte geldt (halve roostercel verschoven)

xbel0 = 0.30;
L  = 29.7;
ib = 4e-4;
C  = 45;
B  = 0.1;
hb = -1*(-6e-5-0.40188);   % waterdiepte op benedenstroomse rand

Q = B * hb * C * sqrt( hb * ib )
N = 100;      % for Belanger-computation

rekendx     = 0.00001;
[xbel,hbel] = belanger2(Q,C,L,B,ib,xb,hb,N,rekendx);
xbel=xbel';
hbel=hbel';
%xbel=xbel+xbel0;

figure;
subplot(2,1,1);
hold on; grid on; box on;
pp=plot(x,d_2Dfull,'-xr');
pp=plot(x,d_2D1cel,'-ob');
pp=plot(mesh1d_node_x,d_1D,'-dg');
pp=plot(xbel,hbel,'-k'); set(pp,'LineWidth',2);
legend('2D','quasi2D','1D','est. uniform','location','best');
xlabel('x [m]')
ylabel('water depth [m]')
xlim([0 31])
title(strrep(caseid,'_','\_'));

subplot(2,1,2);
hold on; grid on; box on;
pp=plot(x,d_2Dfull-hbel(2:end-1),'-xr');
pp=plot(x,d_2D1cel-hbel(2:end-1),'-ob');
pp=plot(mesh1d_node_x,d_1D-hbel(2:end-1),'-dg');
legend('2D','quasi2D','1D','location','best');
xlabel('x [m]')
ylabel('diff. water depth [m]')
xlim([0 31])

x0=11;
y0=0;
width=16.5;
height=16.5;
set(gcf,'units','centimeters','position',[x0,y0,width,height]);
set(gcf,'PaperPositionMode','Auto','PaperUnits','centimeters','PaperSize',[width, height]);  
figname = caseid;
print([figdir,'/pdf/',figname,'.pdf'],'-dpdf','-r300')
print([figdir,figname,'.png'],'-dpng','-r300')

