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

clc;
clear all;
close all;

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

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

nc_sobek3   = ['..\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');
sobek3_bl   = sobek3_wl - sobek3_wd;

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_mor_bl  = ncread(nc_dflowfm,'mesh1d_mor_bl');

%     mesh1d_mor_bl              
%            Size:       41x121
%            Dimensions: nmesh1d_node,time
%            Datatype:   double
%            Attributes:
%                        mesh          = 'mesh1d'
%                        location      = 'node'
%                        coordinates   = 'mesh1d_node_x mesh1d_node_y'
%                        standard_name = ''
%                        long_name     = 'Time-varying bottom level in flow cell center'
%                        units         = 'm'
%                        grid_mapping  = 'projected_coordinate_system'
%                        _FillValue    = -999

%% bovenstrooms al tijdserie

timehr = time'./(3600*24);

figure (1);

subplot(2,1,1);
hold on; box on; grid on;
plot(timehr,sobek3_wl( 2,:),'-xr');
plot(timehr,mesh1d_wl( 2,:),'-ob');
xlabel('time [days]')
ylabel('water level')
ylim([5 6])
title('us')
legend('SOBEK','FM','location','best')

subplot(2,1,2);
hold on; box on; grid on;
plot(timehr,sobek3_wl(end-1,:),'-xr');
plot(timehr,mesh1d_wl(end-1,:),'-ob');
xlabel('time [days]')
ylabel('water level')
ylim([4 5])
title('ds')

figname = [caseid '_fig01_waterlevel_us_ds_Ft'];
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')

%% bed level (x) @ T=0 and T=end

profdata = [ ...
    0	1    ;
4000	0.6  ;
4500	0.95 ;
5000	0.9  ;
5500	0.45 ;
10000	0    ;
];

figure (2);

subplot(1,1,1);
hold on; box on; grid on;
%p=plot( profdata(:,1) + sobek3_x(1) , profdata(:,2)   ,' sg'); set(p,'markerfacecolor','g','markeredgecolor','k','markersize',6);
p=plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,1),'-xr'); set(p,'markersize',8);
p=plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,1),'-ob');
xlabel('x [m]')
ylabel('bed level [m]')
ylim([-0.05 1.05])
title('T = 0 days')
%legend('SOBEK prof locs','SOBEK','FM','location','best')
legend('SOBEK','FM','location','best')

figname = [caseid '_fig02_bedlevel_T0'];
x0=11; y0=0; width=20; height=25;
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')

figure (3);

subplot(1,5,1);
hold on; box on; grid on;
p=plot( profdata(:,1) + sobek3_x(1) , profdata(:,2)   ,' sg'); set(p,'markerfacecolor','g','markeredgecolor','k','markersize',6);
p=plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,1),'-xr'); set(p,'markersize',8);
p=plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,1),'-ob');
xlabel('x [m]')
ylabel('bed level [m]')
ylim([-0.05 1.05])
title('T = 0 days')
legend('SOBEK','FM','location','best')

subplot(1,5,2);
hold on; box on; grid on;
plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,030+1),'-xr');
plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,030+1),'-ob');
xlabel('x [m]')
%ylabel('bed level [m]')
ylim([0 1])
title('T = 30 days')
legend('SOBEK','FM','location','best')

subplot(1,5,3);
hold on; box on; grid on;
plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,060+1),'-xr');
plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,060+1),'-ob');
xlabel('x [m]')
%ylabel('bed level [m]')
ylim([0 1])
title('T = 60 days')
legend('SOBEK','FM','location','best')

subplot(1,5,4);
hold on; box on; grid on;
plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,090+1),'-xr');
plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,090+1),'-ob');
xlabel('x [m]')
%ylabel('bed level [m]')
ylim([0 1])
title('T = 90 days')
legend('SOBEK','FM','location','best')

subplot(1,5,5);
hold on; box on; grid on;
plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,120+1),'-xr');
plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,120+1),'-ob');
xlabel('x [m]')
%ylabel('bed level [m]')
ylim([0 1])
title('T = 120 days')
legend('SOBEK','FM','location','best')
 
figname = [caseid '_fig03_bedlevel_Fxt_per_30_days'];
x0=11; y0=0; width=30; 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')

%% bed level (x,t)
figure (4);

subplot(1,2,1);
hold on; box on; grid on;
for i=1:10:121;
    plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,i) );
end
xlabel('x [m]')
ylabel('bed level [m]')
ylim([0 1])
title('SOBEK3')
%title('T = 0')
%legend('SOBEK','FM','location','best')

subplot(1,2,2);
hold on; box on; grid on;
for i=1:10:121;
    plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,i));
end
%plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,1),'-xb');
ylim([0 1])
xlabel('x [m]')
ylabel('bed level [m]')
title('FM')

figname = [caseid '_fig04_bedlevel_Fx_T0_Tend'];
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')

%% bed level (x,t)

T = 1:20:121;
clr = ['k','r','g','b','m','c','y'];
legstr = {};
pp = [];

figure (5);
subplot(1,1,1);
hold on; box on; grid on;
legt = 1;
for i=1:length(T);
    mapnr = T(i);
       plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,mapnr) , ['--x' clr(i)]);
    p= plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,mapnr) , ['- o' clr(i)]);
    pp = [pp p];
    legstr{i,1} = ['T' num2str(mapnr-1)];
end
xlabel('x [m]')
ylabel('bed level [m]')
ylim([-0.05 1.05])
title('SOBEK3 (--x) vs DFLOWFM (-o)')
legend(pp,legstr)
%title('T = 0')
%legend('SOBEK','FM','location','best')

figname = [caseid '_fig05_bedlevel_Fxt_compare'];
x0=11; y0=0; width=25; height=15;
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')

%% bed level (x,t)

T = 1:30:121;
clr = ['m','r','g','b','k'];

figure (6);

legstr = {};
pp = [];

subplot(2,1,1);
hold on; box on; grid on;
legt = 1;
for i=1:length(T);
    mapnr = T(i);
    p  = plot( sobek3_x( 1:end-0) , sobek3_bl    ( 1:end-0,mapnr) , ['-x' clr(i)]);
    pp = [pp p];
    legstr{i,1} = ['T' num2str(mapnr-1)];
end
xlabel('x [m]')
ylabel('bed level [m]')
ylim([-0.05 1.05])
title('SOBEK3' )
legend(pp,legstr)

legstr = {};
pp = [];

subplot(2,1,2);
hold on; box on; grid on;
legt = 1;
for i=1:length(T);
    mapnr = T(i);
    p  = plot( mesh1d_x( 1:end-0) , mesh1d_mor_bl( 1:end-0,mapnr) , ['-x' clr(i)]);
    pp = [pp p];
    legstr{i,1} = ['T' num2str(mapnr-1)];
end
xlabel('x [m]')
ylabel('bed level [m]')
ylim([-0.05 1.05])
title('DFLOWFM' )
legend(pp,legstr)
figname = [caseid '_fig06_bedlevel_Fxt_compare'];
x0=11; y0=0; width=20; height=20;
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')

