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

clc;
clear all;
close all;

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

plot_hydr = true;  %make figures

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

nc_sobek3   = ['..\sobek3_dimr\dflow1d\output\gridpoints.nc'];
nc_dflowfm  = ['..\dflowfmoutput\' caseid '_map.nc'];

info=ncinfo(nc_sobek3);
ncdisp(nc_sobek3);
time=ncread(nc_sobek3,'time');

sobek3_x    = ncread(nc_sobek3,'x_coordinate');
sobek3_y    = ncread(nc_sobek3,'y_coordinate');
sobek3_wd   = ncread(nc_sobek3,'water_depth');
sobek3_wl   = ncread(nc_sobek3,'water_level');

info=ncinfo(nc_dflowfm);
ncdisp(nc_dflowfm);
time=ncread(nc_dflowfm,'time');

mesh1d_node_x  = ncread(nc_dflowfm,'mesh1d_node_x');
mesh1d_node_y  = ncread(nc_dflowfm,'mesh1d_node_y');
mesh1d_x       = mesh1d_node_x;
mesh1d_y       = mesh1d_node_y;
mesh1d_wd      = ncread(nc_dflowfm,'mesh1d_waterdepth');
mesh1d_wl      = ncread(nc_dflowfm,'mesh1d_s1');
mesh1d_d50     = ncread(nc_dflowfm,'mesh1d_DXX04');
mesh1d_d90     = ncread(nc_dflowfm,'mesh1d_DXX06');

mesh1d_sxtot   = squeeze(ncread(nc_dflowfm,'mesh1d_sxtot'));
mesh1d_sytot   = squeeze(ncread(nc_dflowfm,'mesh1d_sytot'));

% let op: sbcx is de juiste, nog niet getransformeerde transport waar
% we in dit geval mee moeten vergelijken.
mesh1d_sbcx    = squeeze(ncread(nc_dflowfm,'mesh1d_sbcx'));

mesh1d_uuu     = squeeze(ncread(nc_dflowfm,'mesh1d_uuu'));

% parameters needed to compute sediment transport (kg/s)
mesh1d_czs = squeeze(ncread(nc_dflowfm,'mesh1d_czs')); 
C = mesh1d_czs(:,1);
mesh1d_d50   = squeeze(ncread(nc_dflowfm,'mesh1d_DXX04'));
D50 = mesh1d_d50(:,1);

B = 200;
g = 9.81;
rhow = 1000;
rhos = 2650;
delt = (rhos-rhow)/rhow;
alph = 1;

%% bovenstrooms als tijdserie

if plot_hydr == 1;

    timehr = time'./3600;

    figure;
    subplot(2,2,1);
    hold on; box on; grid on;
    plot(timehr,sobek3_wl( 2,:),'-xr');
    plot(timehr,mesh1d_wl( 2,:),'-ob');
    %xlabel('time [hr]')
    ylabel('water level')
    title('T1')
    legend('SOBEK','FM','location','best')

    subplot(2,2,2);
    hold on; box on; grid on;
    plot(timehr,sobek3_wl(12,:),'-xr');
    plot(timehr,mesh1d_wl(12,:),'-ob');
    %xlabel('time [hr]')
    %ylabel('water level')
    title('T2')

    subplot(2,2,3);
    hold on; box on; grid on;
    plot(timehr,sobek3_wl(22,:),'-xr');
    plot(timehr,mesh1d_wl(22,:),'-ob');
    xlabel('time [hr]')
    ylabel('water level')
    title('T3')

    subplot(2,2,4);
    hold on; box on; grid on;
    plot(timehr,sobek3_wl(32,:),'-xr');
    plot(timehr,mesh1d_wl(32,:),'-ob');
    xlabel('time [hr]')
    %ylabel('water level')
    title('T4')

    figname = [caseid '_fig01_waterlevelnode2T1'];
    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]);
    %print([figdir,'/pdf/',figname,'.pdf'],'-dpdf','-r300')
    print([figdir,figname,'.png'],'-dpng','-r300')

    %% verschil T1

    figure;
    hold on; box on; grid on;
    plot(timehr,mesh1d_wl( 2,:) - sobek3_wl( 2,:),'-xb');
    xlabel('time [hr]')
    ylabel('water level')
    title('T1')
    legend('FM minus SOBEK','location','best')

    figname = [caseid '_fig02_waterleveldifferencenode2T1'];
    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]);
    %print([figdir,'/pdf/',figname,'.pdf'],'-dpdf','-r300')
    print([figdir,figname,'.png'],'-dpng','-r300')

    %% water level(x) T1

    figure;
    subplot(3,1,1);
    hold on; box on; grid on;
    plot(sobek3_x( 2:10),sobek3_wl( 2:10,24+1),'-xr');
    plot(mesh1d_x( 1:10),mesh1d_wl( 1:10,24+1),'-xb');
    xlabel('x [m]')
    ylabel('water level')
    title('T1 - after 1 day')
    legend('SOBEK','FM','location','best')

    subplot(3,1,2);
    hold on; box on; grid on;
    plot(sobek3_x( 2:10),sobek3_wl( 2:10,48+1),'-xr');
    plot(mesh1d_x( 1:10),mesh1d_wl( 1:10,48+1),'-xb');
    xlabel('x [m]')
    ylabel('water level')
    title('T1 - after 2 days')
    legend('SOBEK','FM','location','best')

    subplot(3,1,3);
    hold on; box on; grid on;
    plot(sobek3_x( 2:10),sobek3_wl( 2:10,72+1),'-xr');
    plot(mesh1d_x( 1:10),mesh1d_wl( 1:10,72+1),'-xb');
    xlabel('x [m]')
    ylabel('water level')
    title('T1 - after 3 days')
    legend('SOBEK','FM','location','best')

    figname = [caseid '_fig03_waterlevelT1'];
    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]);
    %print([figdir,'/pdf/',figname,'.pdf'],'-dpdf','-r300')
    print([figdir,figname,'.png'],'-dpng','-r300')

    %% water level difference (x) T1

    figure;
    subplot(3,1,1);
    hold on; box on; grid on;
    plot(mesh1d_x( 2:10),mesh1d_wl( 2:10,24+1) - sobek3_wl( 2:10,24+1),'-xb');
    xlabel('x [m]')
    ylabel('water level diff [m]')
    title('T1 - after 1 day')
    legend('FM minus SOBEK','location','best')

    subplot(3,1,2);
    hold on; box on; grid on;
    plot(mesh1d_x( 2:10),mesh1d_wl( 2:10,48+1) - sobek3_wl( 2:10,48+1),'-xb');
    xlabel('x [m]')
    ylabel('water level diff [m]')
    title('T1 - after 2 days')
    legend('FM minus SOBEK','location','best')

    subplot(3,1,3);
    hold on; box on; grid on;
    plot(mesh1d_x( 2:10),mesh1d_wl( 2:10,72+1) - sobek3_wl( 2:10,72+1),'-xb');
    xlabel('x [m]')
    ylabel('water level diff [m]')
    title('T1 - after 3 days')
    legend('FM minus SOBEK','location','best')

    figname = [caseid '_fig04_waterleveldifferenceT1'];
    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]);
    %print([figdir,'/pdf/',figname,'.pdf'],'-dpdf','-r300')
    print([figdir,figname,'.png'],'-dpng','-r300')

    %% water depth(x) T1

    figure;
    subplot(3,1,1);
    hold on; box on; grid on;
    plot(sobek3_x( 1:10),sobek3_wd( 1:10,24+1),'-xr');
    plot(mesh1d_x( 1:10),mesh1d_wd( 1:10,24+1),'-xb');
    xlabel('x [m]')
    ylabel('water depth')
    title('T1 - after 1 day')
    legend('SOBEK','FM','location','best')

    subplot(3,1,2);
    hold on; box on; grid on;
    plot(sobek3_x( 1:10),sobek3_wd( 1:10,48+1),'-xr');
    plot(mesh1d_x( 1:10),mesh1d_wd( 1:10,48+1),'-xb');
    xlabel('x [m]')
    ylabel('water depth')
    title('T1 - after 2 days')
    legend('SOBEK','FM','location','best')

    subplot(3,1,3);
    hold on; box on; grid on;
    plot(sobek3_x( 1:10),sobek3_wd( 1:10,72+1),'-xr');
    plot(mesh1d_x( 1:10),mesh1d_wd( 1:10,72+1),'-xb');
    xlabel('x [m]')
    ylabel('water depth')
    title('T1 - after 3 days')
    legend('SOBEK','FM','location','best')

    figname = [caseid '_fig05_waterdepthT1'];
    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]);
    %print([figdir,'/pdf/',figname,'.pdf'],'-dpdf','-r300')
    print([figdir,figname,'.png'],'-dpng','-r300')

end

%% sediment transport

T1 = 24;
T2 = 48;
T3 = 72;

%Ureken = mesh1d_ucxq;
Ureken = mesh1d_uuu;

% parameters
Acal = 1.0;
shld_cr = 0.047;

D90 = mesh1d_d90( :, 1);
D50 = mesh1d_d50( :, 1);

u_Q1_fm = Ureken( :, T1+1);
h_Q1_fm = mesh1d_wd( :, T1+1);
% calculate ripple factor
C90 = 18*log10(12*h_Q1_fm./D90);
mu = min((C./C90).^1.5,1);
xi = 1; %hiding exposure coefficient
shld = (u_Q1_fm./C).^2./delt./D50;

S_Q1_analytic = B.*rhos .* meyer_peter_muller(Acal,D50,delt,g,mu,shld,xi,shld_cr);  % kg/s
S_Q1_fm       = mesh1d_sbcx( :, T1+1) .* rhos .* B;                                % kg/s

u_Q2_fm = Ureken( :, T2+1);
h_Q2_fm = mesh1d_wd( :, T2+1);
% calculate ripple factor
C90 = 18*log10(12*h_Q2_fm./D90);
mu = min((C./C90).^1.5,1);
xi = 1; %hiding exposure coefficient
shld = (u_Q2_fm./C).^2./delt./D50;

S_Q2_analytic = B.*rhos .* meyer_peter_muller(Acal,D50,delt,g,mu,shld,xi,shld_cr);  % kg/s
S_Q2_fm       = mesh1d_sbcx( :, T2+1) .* rhos .* B;                                % kg/s

u_Q3_fm       = Ureken( :, T3+1);
h_Q3_fm = mesh1d_wd( :, T3+1);
% calculate ripple factor
C90 = 18*log10(12*h_Q3_fm./D90);
mu = min((C./C90).^1.5,1);
xi = 1; %hiding exposure coefficient
shld = (u_Q3_fm./C).^2./delt./D50;

S_Q3_analytic = B.*rhos .* meyer_peter_muller(Acal,D50,delt,g,mu,shld,xi,shld_cr);  % kg/s
S_Q3_fm       = mesh1d_sbcx( :, T3+1) .* rhos .* B;                                % kg/s

%% sediment transport naar tabel schrijven

xlsfile = [figdir 'sedtr.xlsx'];

%return

%xlswrite(FILE,ARRAY,SHEET,RANGE)
%
xlswrite (xlsfile,S_Q1_fm      ( 1:10)','Sheet1','C4');
xlswrite (xlsfile,S_Q1_fm      (11:20)','Sheet1','C5');
xlswrite (xlsfile,S_Q1_fm      (21:30)','Sheet1','C6');
xlswrite (xlsfile,S_Q1_fm      (31:40)','Sheet1','C7');
%
xlswrite (xlsfile,S_Q1_analytic( 1:10)','Sheet1','C10');
xlswrite (xlsfile,S_Q1_analytic(11:20)','Sheet1','C11');
xlswrite (xlsfile,S_Q1_analytic(21:30)','Sheet1','C12');
xlswrite (xlsfile,S_Q1_analytic(31:40)','Sheet1','C13');
%
xlswrite (xlsfile,S_Q2_fm      ( 1:10)','Sheet1','O4');
xlswrite (xlsfile,S_Q2_fm      (11:20)','Sheet1','O5');
xlswrite (xlsfile,S_Q2_fm      (21:30)','Sheet1','O6');
xlswrite (xlsfile,S_Q2_fm      (31:40)','Sheet1','O7');
%
xlswrite (xlsfile,S_Q2_analytic( 1:10)','Sheet1','O10');
xlswrite (xlsfile,S_Q2_analytic(11:20)','Sheet1','O11');
xlswrite (xlsfile,S_Q2_analytic(21:30)','Sheet1','O12');
xlswrite (xlsfile,S_Q2_analytic(31:40)','Sheet1','O13');
%
xlswrite (xlsfile,S_Q3_fm      ( 1:10)','Sheet1','AA4');
xlswrite (xlsfile,S_Q3_fm      (11:20)','Sheet1','AA5');
xlswrite (xlsfile,S_Q3_fm      (21:30)','Sheet1','AA6');
xlswrite (xlsfile,S_Q3_fm      (31:40)','Sheet1','AA7');
%
xlswrite (xlsfile,S_Q3_analytic( 1:10)','Sheet1','AA10');
xlswrite (xlsfile,S_Q3_analytic(11:20)','Sheet1','AA11');
xlswrite (xlsfile,S_Q3_analytic(21:30)','Sheet1','AA12');
xlswrite (xlsfile,S_Q3_analytic(31:40)','Sheet1','AA13');
%
