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

clc;
clear all;
close all;

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

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

simdirs = {'c20_straight_channel_1Dvs2D_zero_bedslope/',...
           'c21_straight_channel_1Dvs2D_zero_bedslope_blev3/',...
           'c22_straight_channel_1Dvs2D_zero_bedslope_blev6/'};
       
caseids = {'t08_belangerflat_bedlevtype1', 't08_belangerflat_bedlevtype3', 't08_belangerflat_bedlevtype6'}; 
       
for k = 1:length(caseids); 
    caseid  = caseids{k};

    datadir = ['../../',simdirs{k},'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);

    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);

    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= caseid;
    id_time = -1;
    P03_plotoption3;

    %% manipulate water level data

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

    % selecteren centerline uit 2D grid
    id_F2D    = find( mesh2d_face_y>1.30 & mesh2d_face_y<1.40 );
    mesh2d_face_x      = mesh2d_face_x(id_F2D);
    mesh2d_face_y      = mesh2d_face_y(id_F2D);
    mesh2d_waterdepth  = mesh2d_waterdepth(id_F2D);
    clear id_F2D;

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

    d_1D     = mesh1d_waterdepth;
    d_2D     = mesh2d_waterdepth;

    %% belanger

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

    ib = 0;
    C  = 45;
    B  = 0.3;
    L  = 30.;

    %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);
    rekendx     = 0.00001;
    [xbel,hbel] = belanger2(Q,C,L,B,ib,xb,hb,N,rekendx);
    xbel=xbel';
    hbel=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','1D','analytic','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(mesh2d_face_x,d_2D-hbel(2:end-1),'-xr');
    pp=plot(mesh1d_node_x,d_1D-hbel(2:end-1),'-dg');
    legend('2D','1D','location','best');
    xlabel('x [m]')
    ylabel('water depth diff wrt analytic [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')
end
