%% 2D Dam Break Plot Generator
% Tim Leijnse
% Made to generate plots for the 2D Dam Break test as performed by Stelling and Duinmeijer (2003)

function plot_2Ddrybed(fname, dataDir, dx, dy, nmax, zs_sb, zst_sb, zs_nh, zst_nh, dtout, timesteps)

%% figure h3    
fh3 = figure('units','pixels','position',[0 0 1920 1080]); 

%Actions for Stelling picture
C = imread(fullfile( dataDir, 'Stelling_n0012_Dry.png'));
B = imcrop(C,[253 27 1016-253 430-27]);
x = (1:1:764) -68;
y = fliplr((1:1:404)-53);
[X,Y] = meshgrid(x,y);
Z = X;
Z(:,:) = 0;
factor_x = 686/8 ;
factor_y = 342/4 ;

surf(X,Y,Z,B,'EdgeColor', 'none');
view(2)
grid off
axis off
hold on

%For scatterplot
t = [1+1,2+1,3+1,4+1]; %time steps you want to show
xxx = 0:dy*factor_y:8*factor_y;

index_first = zeros(length(t),nmax+1);
index_last = zeros(length(t),nmax+1);

for kk=1:length(t)
    zs_front_temp = squeeze( zs_sb(t(kk),:,26:end)>0);  % % location 26 depends on grid size

    for ll=1:nmax+1
        last = find(zs_front_temp(ll,:),1,'last');
        first = find(zs_front_temp(ll,:),1,'first');
        if isempty(last)
            last = NaN;
        end
        if isempty(first)
            first = NaN;
        end
        index_last(kk,ll) = last;
        index_first(kk,ll) = first;

    end
    scatter((index_last(kk,41:end)*dx*factor_x)-dx*factor_x,xxx(41:end)-4*factor_y,'ro','filled','DisplayName','zs_1'); 
    hold on
    scatter((index_first(kk,41:end)*dx*factor_x)-dx*factor_x,xxx(41:end)-4*factor_y,'ro','filled','DisplayName','zs_1'); 
    hold on
end

index_first = zeros(length(t),nmax+1);
index_last = zeros(length(t),nmax+1);

for kk=1:length(t)
    zs_front_temp = squeeze( zs_nh(t(kk),:,26:end)>0);  % % location 26 depends on grid size

    for ll=1:nmax+1
        last = find(zs_front_temp(ll,:),1,'last');
        first = find(zs_front_temp(ll,:),1,'first');
        if isempty(last)
            last = NaN;
        end
        if isempty(first)
            first = NaN;
        end
        index_last(kk,ll) = last;
        index_first(kk,ll) = first;

    end
    scatter((index_last(kk,41:end)*dx*factor_x)-dx*factor_x,xxx(41:end)-4*factor_y,'bo','filled','DisplayName','zs_1'); 
    hold on
    scatter((index_first(kk,41:end)*dx*factor_x)-dx*factor_x,xxx(41:end)-4*factor_y,'bo','filled','DisplayName','zs_1'); 
    hold on
end

hold off

xb_write_plot(fh3, fname, '03_DamBreak2Ddrybed');


%% figure h4

zst_time = 0:(dtout/10):(timesteps-1);
zst_sb_loc1_5 = zst_sb(:,2-1)-0.3;
zst_sb_loc3_5 = zst_sb(:,3-1);
zst_sb_loc8_5 = zst_sb(:,4-1);
zst_sb_loc11_5 = zst_sb(:,5-1);

zst_nh_loc1_5 = zst_nh(:,2-1)-0.3;
zst_nh_loc3_5 = zst_nh(:,3-1);
zst_nh_loc8_5 = zst_nh(:,4-1);
zst_nh_loc11_5 = zst_nh(:,5-1);        

%Limit values to axis limits
zst_sb_loc1_5(zst_sb_loc1_5 > 0.3) = 0.3;
zst_sb_loc1_5(zst_sb_loc1_5 < 0) = 0;
zst_sb_loc3_5(zst_sb_loc3_5 > 0.08) = 0.08;
zst_sb_loc3_5(zst_sb_loc3_5 < 0) = 0;
zst_sb_loc8_5(zst_sb_loc8_5 > 0.06) = 0.06;
zst_sb_loc8_5(zst_sb_loc8_5 < 0) = 0;
zst_sb_loc11_5(zst_sb_loc11_5 > 0.06) = 0.06;
zst_sb_loc11_5(zst_sb_loc11_5 < 0) = 0;

zst_nh_loc1_5(zst_nh_loc1_5 > 0.3) = 0.3;
zst_nh_loc1_5(zst_nh_loc1_5 < 0) = 0;
zst_nh_loc3_5(zst_nh_loc3_5 > 0.08) = 0.08;
zst_nh_loc3_5(zst_nh_loc3_5 < 0) = 0;
zst_nh_loc8_5(zst_nh_loc8_5 > 0.06) = 0.06;
zst_nh_loc8_5(zst_nh_loc8_5 < 0) = 0;
zst_nh_loc11_5(zst_nh_loc11_5 > 0.06) = 0.06;
zst_nh_loc11_5(zst_nh_loc11_5 < 0) = 0;

fh4 = figure('units','pixels','position',[0 0 1920 1080]); 
%Actions for Stelling picture
%Gauge 00
C = imread(fullfile(dataDir, 'Stelling_Dry_Gauge00and01.png'));
B0 = imcrop(C,[176 34 611-176 406-34]);
x0 = (1:1:436) -49;
y0 = fliplr((1:1:373)-87);
[X0,Y0] = meshgrid(x0,y0);
Z0 = X0;
Z0(:,:) = 0;
factor_x0 = 382/26;
factor_y0 = 276/0.3;
%Gauge 01
B1 = imcrop(C,[642 34 1085-642 406-34]);
x1 = (1:1:444) -53;
y1 = fliplr((1:1:373)-87);
[X1,Y1] = meshgrid(x1,y1);
Z1 = X1;
Z1(:,:) = 0;
factor_x1 = 382/26;
factor_y1 = 276/0.08;
%Gauge 02
C2 = imread(fullfile(dataDir, 'Stelling_Dry_Gauge02and03.png'));
B2 = imcrop(C2,[185 22 626-185 385-22]);
x2 = (1:1:442) -53;
y2 = fliplr((1:1:364)-80);
[X2,Y2] = meshgrid(x2,y2);
Z2 = X2;
Z2(:,:) = 0;
factor_x2 = 382/26;
factor_y2 = 276/0.06;
%Gauge 03
B3 = imcrop(C2,[663 22 1104-663 385-22]);
x3 = (1:1:442) -48;
y3 = fliplr((1:1:364)-80);
[X3,Y3] = meshgrid(x3,y3);
Z3 = X3;
Z3(:,:) = 0;
factor_x3 = 382/26;
factor_y3 = 276/0.06;

%The actual plots
subplot('Position',[0.05 0.5 0.45 0.45])
surf(X0,Y0,Z0,B0,'EdgeColor', 'none');
view(2)
grid off
axis off
hold on
plot(zst_time*factor_x0,zst_sb_loc1_5*factor_y0,'k','LineWidth',4) 
hold on
plot(zst_time*factor_x0,zst_nh_loc1_5*factor_y0,'cyan','LineWidth',4) 
hold on

subplot('Position',[0.5 0.5 0.5 0.45])
surf(X1,Y1,Z1,B1,'EdgeColor', 'none');
view(2)
grid off
axis off
hold on
plot(zst_time*factor_x1,zst_sb_loc3_5*factor_y1,'k','LineWidth',4)
hold on
plot(zst_time*factor_x1,zst_nh_loc3_5*factor_y1,'cyan','LineWidth',4)
hold on

subplot('Position',[0 0.05 0.5 0.45])
surf(X2,Y2,Z2,B2,'EdgeColor', 'none');
view(2)
grid off
axis off
hold on    
plot(zst_time*factor_x2,zst_sb_loc8_5*factor_y2,'k','LineWidth',4)
hold on
plot(zst_time*factor_x2,zst_nh_loc8_5*factor_y2,'cyan','LineWidth',4)
hold on

subplot('Position',[0.55 0.05 0.45 0.45])
surf(X3,Y3,Z3,B3,'EdgeColor', 'none');
view(2)
grid off
axis off
hold on
plot(zst_time*factor_x3,zst_sb_loc11_5*factor_y3,'k','LineWidth',4)
hold on
plot(zst_time*factor_x3,zst_nh_loc11_5*factor_y3,'cyan','LineWidth',4)

hold off
xb_write_plot(fh4, fname, '04_DamBreak2Ddrybed');