% 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

simdirs = {'c30_Channels_1D_zigzag_zero_bedslope/', ...
           'c31_Channels_1D_zigzag_zero_bedslope_blev3/'};

caseids = {'t07_conv2Dmin1_blev1',...
           't07_conv2Dmin1_blev3'};

for k = 1:length(caseids);
       
    figdir  = '../doc/figures/';
    caseid  = caseids{k};

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

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

    % selecteren rechte lijn 1D grid
    id_1D_s     = find( mesh1d_node_y<3.50 );
    mesh1d_node_x_s      = mesh1d_node_x(id_1D_s);
    mesh1d_node_y_s      = mesh1d_node_y(id_1D_s);
    mesh1d_waterdepth_s  = mesh1d_waterdepth(id_1D_s);
    clear id_1D_s;

    [mesh1d_node_x_s,I] = sortrows(mesh1d_node_x_s);
    mesh1d_node_y_s     = mesh1d_node_y_s(I);
    mesh1d_waterdepth_s = mesh1d_waterdepth_s(I);
    clear I;

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

    [mesh2d_face_x_s,I] = sortrows(mesh2d_face_x_s);
    mesh2d_waterdepth_s = mesh2d_waterdepth_s(I);
    clear I;

    d_1D     = mesh1d_waterdepth_s;
    d_2D     = mesh2d_waterdepth_s;
    x_1D     = mesh1d_node_x_s;
    x_2D     = mesh2d_face_x_s;

    %% 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(x_2D,d_2D,'-xr');
    pp=plot(x_1D,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(x_2D,d_2D-hbel(2:end-1),'-xr');
    pp=plot(x_1D,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')

    %% nu ook voor de channels met bends

    % 2D with single bend
    id_2Db              = find( mesh2d_face_y>4.00 & mesh2d_face_x<15.5 );
    mesh2d_face_x_b     = mesh2d_face_x(id_2Db);
    mesh2d_face_y_b     = mesh2d_face_y(id_2Db);
    mesh2d_waterdepth_b = mesh2d_waterdepth(id_2Db);
    clear id_2Db;

    % sort based on x AND y, first X, then ascend Y
    xy = [mesh2d_face_x_b mesh2d_face_y_b];
    [xy,I] = sortrows(xy);
    mesh2d_face_x_b     = mesh2d_face_x_b(I);
    mesh2d_waterdepth_b = mesh2d_waterdepth_b(I);
    clear xy I;

    % 1D with single bend
    id_1D_b     = find( mesh1d_node_y>3.50 & mesh1d_node_x>15.5 );
    mesh1d_node_x_b      = mesh1d_node_x(id_1D_b);
    mesh1d_node_y_b      = mesh1d_node_y(id_1D_b);
    mesh1d_waterdepth_b  = mesh1d_waterdepth(id_1D_b);
    clear id_1D_b;

    xy = [mesh1d_node_x_b mesh1d_node_y_b];
    xy(:,1) = round(xy(:,1),2);
    [xy,I] = sortrows(xy,[1 2]);
    mesh1d_node_x_b     = mesh1d_node_x_b(I);
    mesh1d_node_y_b     = mesh1d_node_y_b(I);
    mesh1d_waterdepth_b = mesh1d_waterdepth_b(I);
    clear xy I;

    % 1D with zigzag
    id_1D_z     = find( mesh1d_node_y>3.50 & mesh1d_node_x<15.5 );
    mesh1d_node_x_z      = mesh1d_node_x(id_1D_z);
    mesh1d_node_y_z      = mesh1d_node_y(id_1D_z);
    mesh1d_waterdepth_z  = mesh1d_waterdepth(id_1D_z);
    clear id_1D_z;

    xy = [mesh1d_node_x_z mesh1d_node_y_z];
    xy = round(xy,3);
    [xy,I] = sortrows(xy,[1 2]);
    [mesh1d_node_x_z,I] = sortrows(mesh1d_node_x_z);
    mesh1d_node_y_z     = mesh1d_node_y_z(I);
    mesh1d_waterdepth_z = mesh1d_waterdepth_z(I);
    clear I;

    d_2D_b = mesh2d_waterdepth_b;
    d_1D_b = mesh1d_waterdepth_b;
    d_1D_z = mesh1d_waterdepth_z;

    figure;
    subplot(1,1,1);
    hold on; grid on; box on;
    pp=plot(xbel,hbel,'- k'); %set(pp,'LineWidth',2);
    pp=plot(x_2D,d_2D_b,'-xr');
    pp=plot(x_1D,d_1D_b,'-dg');
    pp=plot(x_1D,d_1D_z,'-vb');
    legend('semi-analytic','2D - 1 bend','1D - 1 bend','1D - multiple bends','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_2D,d_2D-hbel(2:end-1),'-xr');
    % pp=plot(x_1D,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, '_bends'];
    %print([figdir,'/pdf/',figname,'.pdf'],'-dpdf','-r300')
    print([figdir,figname,'.png'],'-dpng','-r300')

end