%% ART4_Plants_to_Algae_bifurc_FIG105.m

% Roeland C. van de Vijsel
% 27/05/2020
% NIOZ Yerseke, The Netherlands
% www.nioz.nl

% This Matlab script analyses creates plots for Figures 5.6 and 5.7 in Chapter 5 of the
% PhD Thesis of Roeland C. van de Vijsel (https://doi.org/10.33612/diss.160081233).

% The only thing that needs to be manually varied in this script, is number 
% assigned to "SUBRUN2" (line 31): choose from numbers 1,2,4,11,12 and 14.

clear all
close all
clc

%%
MORF = 44712;

plotting = 0;
el_end = 41;
logplot = 0;
shade = 1;
slope = 0.25/200; %m/m
nb = 0.016;
nv = 0.2;
kk = 1500;
ampvarmea = 4; % 3 = abs with min/max range; 4 = std
set(0,'defaultfigurecolor',[1 1 1])
for SUBRUN2 = 14%1:25;

for FIG = [105]%41
    orient('landscape')
    close(figure(101))        

        addpath(['FIG',num2str(FIG)])

        Hin_0 = 1.4e-5;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if (FIG == 2) 
            Hin_arr = [1:10,9:-1:1]*Hin_0;
            xmin = min(Hin_arr); xmax = max(Hin_arr)/2;
            D0 = 1e-7*MORF*1;
            EndTime = 15000;
            NX = 512;
            NY = NX;
            pert = 'perturbB';
        elseif (FIG == 3)
            Hin_arr = [1:10,9:-1:1]*Hin_0;
            xmin = min(Hin_arr); xmax = max(Hin_arr)/2;
            D0 = 1e-7*MORF*100;
            EndTime = 15000;
            NX = 512;
            NY = NX;
            pert = 'perturbB';
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        elseif (FIG == 4) 
            Hin_arr = [1:0.5:10,10:-0.5:1,1/2,1/4,1/8,1/16,1/16,1/8,1/4,1/2,1]*Hin_0;
            xmin = min(Hin_arr); xmax = max(Hin_arr)/2;
            Hin_arr = Hin_arr(1:el_end);
            dir_arr = zeros(length(Hin_arr),1); dir_arr;
            D0 = 1e-7*MORF*1;
            EndTime = 15000;
            NX = 512;
            NY = NX;
            el0 = 43; % First element of upward branch
            el1 = 9;  % Last element of upward branch
            pert = 'perturbB';
        elseif (FIG == 5) 
            Hin_arr = [1:0.5:10,10:-0.5:1,1/2,1/4,1/8,1/16,1/16,1/8,1/4,1/2,1]*Hin_0;
            xmin = min(Hin_arr); xmax = max(Hin_arr)/2;
            Hin_arr = Hin_arr(1:el_end);
            D0 = 1e-7*MORF*10;
            EndTime = 15000;
            NX = 512;
            NY = NX;
            el0 = 43; % First element of upward branch
            el1 = 9;  % Last element of upward branch
            pert = 'perturbB';
        elseif (FIG == 6) 
            Hin_arr = [1:0.5:10,10:-0.5:1,1/2,1/4,1/8,1/16,1/16,1/8,1/4,1/2,1]*Hin_0;
            xmin = min(Hin_arr); xmax = max(Hin_arr)/2;
            Hin_arr = Hin_arr(1:el_end);
            D0 = 1e-7*MORF*100;
            EndTime = 15000;
            NX = 512;
            NY = NX;
            el0 = 43; % First element of upward branch
            el1 = 9;  % Last element of upward branch
            pert = 'perturbB';
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        elseif (FIG == 7)
            Hin_arr = [1:0.5:10,10:-0.5:1,1/2,1/4,1/8,1/16,1/16,1/8,1/4,1/2,1]*Hin_0;
            xmin = min(Hin_arr); xmax = max(Hin_arr)/2;
            Hin_arr = Hin_arr(1:el_end);
            dir_arr = zeros(length(Hin_arr),1); dir_arr;
            D0 = 1e-7*MORF*1;
            EndTime = 15000;
            NX = 512;
            NY = NX;
            el0 = 43; % First element of upward branch
            el1 = 9;  % Last element of upward branch
            pert = 'perturbS';
        elseif (FIG == 8)
            Hin_arr = [1:0.5:10,10:-0.5:1,1/2,1/4,1/8,1/16,1/16,1/8,1/4,1/2,1]*Hin_0;
            xmin = min(Hin_arr); xmax = max(Hin_arr)/2;
            Hin_arr = Hin_arr(1:el_end);
            D0 = 1e-7*MORF*10;
            EndTime = 15000;
            NX = 512;
            NY = NX;
            el0 = 43; % First element of upward branch
            el1 = 9;  % Last element of upward branch
            pert = 'perturbS';
        elseif (FIG == 9)
            Hin_arr = [1:0.5:10,10:-0.5:1,1/2,1/4,1/8,1/16,1/16,1/8,1/4,1/2,1]*Hin_0;
            xmin = min(Hin_arr); xmax = max(Hin_arr)/2;
            Hin_arr = Hin_arr(1:el_end);
            D0 = 1e-7*MORF*100;
            EndTime = 15000;
            NX = 512;
            NY = NX;
            el0 = 43; % First element of upward branch
            el1 = 9;  % Last element of upward branch
            pert = 'perturbS';
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         elseif (FIG == 10)
            Hin_arr = [1.0:0.3:5.8, 5.8:-0.3:0.1, 0.1:0.3:1.0]*Hin_0;
            
            Hin_arr = Hin_arr(1:el_end);
            D0 = 1e-7*MORF*1;
            EndTime = 15000;
            NX = 256;
            NY = NX;
            el0 = 38; % First element of upward branch
            el1 = 17;  % Last element of upward branch
            LengthX = 50; % m
            DZ_TF = LengthX*slope*0.5; 
            DZ_MHWS = (MORF*Hin_arr) - DZ_TF;
            Hin_arr = -DZ_MHWS;
            xmin = min(Hin_arr); xmax = max(Hin_arr);
            pert = 'perturbB';
        elseif (FIG == 11)
            Hin_arr = [1.0:0.3:5.8, 5.8:-0.3:0.1, 0.1:0.3:1.0]*Hin_0;
            Hin_arr = Hin_arr(1:el_end);
            D0 = 1e-7*MORF*100;
            EndTime = 15000;
            NX = 256;
            NY = NX;
            el0 = 38; % First element of upward branch
            el1 = 17;  % Last element of upward branch
            LengthX = 50; % m
            DZ_TF = LengthX*slope*0.5; 
            DZ_MHWS = (MORF*Hin_arr) - DZ_TF;
            Hin_arr = -DZ_MHWS;
            xmin = min(Hin_arr); xmax = max(Hin_arr);

            pert = 'perturbB';
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        elseif (FIG == 12)
            Hin_arr = [1.0:0.3:5.8, 5.8:-0.3:0.1, 0.1:0.3:1.0]*Hin_0;
            
            Hin_arr = Hin_arr(1:el_end);
            D0 = 1e-7*MORF*1;
            EndTime = 30000;
            NX = 256;
            NY = NX;
            el0 = 38; % First element of upward branch
            el1 = 17;  % Last element of upward branch
            LengthX = 50; % m
            DZ_TF = LengthX*slope*0.5; 
            DZ_MHWS = (MORF*Hin_arr) - DZ_TF;
            Hin_arr = -DZ_MHWS;
            xmin = min(Hin_arr); xmax = max(Hin_arr);
            pert = 'perturbB';
        elseif (FIG == 13)
            Hin_arr = [1.0:0.3:5.8, 5.8:-0.3:0.1, 0.1:0.3:1.0]*Hin_0;
            
            Hin_arr = Hin_arr(1:el_end);
            D0 = 1e-7*MORF*100;
            EndTime = 30000;
            NX = 256;
            NY = NX;
            el0 = 38; % First element of upward branch
            el1 = 17;  % Last element of upward branch
            LengthX = 50; % m
            DZ_TF = LengthX*slope*0.5; 
            DZ_MHWS = (MORF*Hin_arr) - DZ_TF;
            Hin_arr = -DZ_MHWS;
            xmin = min(Hin_arr); xmax = max(Hin_arr);
            pert = 'perturbB';
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         elseif (FIG == 14)
            Hin_arr = [1.0:0.3:5.8, 5.8:-0.3:0.1, 0.1:0.3:1.0]*Hin_0;
            
            Hin_arr = Hin_arr(1:el_end);
            D0 = 1e-7*MORF*1;
            EndTime = 15000;
            NX = 256;
            NY = NX;
            el0 = 38; % First element of upward branch
            el1 = 17;  % Last element of upward branch
            LengthX = 50; % m
            DZ_TF = LengthX*slope*0.5; 
            DZ_MHWS = (MORF*Hin_arr) - DZ_TF;
            Hin_arr = -DZ_MHWS;
            xmin = min(Hin_arr); xmax = max(Hin_arr);
            pert = 'perturbS';
        elseif (FIG == 15)
            Hin_arr = [1.0:0.3:5.8, 5.8:-0.3:0.1, 0.1:0.3:1.0]*Hin_0;
            
            Hin_arr = Hin_arr(1:el_end);
            D0 = 1e-7*MORF*100;
            EndTime = 15000;
            NX = 256;
            NY = NX;
            el0 = 38; % First element of upward branch
            el1 = 17;  % Last element of upward branch
            LengthX = 50; % m
            DZ_TF = LengthX*slope*0.5; 
            DZ_MHWS = (MORF*Hin_arr) - DZ_TF;
            Hin_arr = -DZ_MHWS;
            xmin = min(Hin_arr); xmax = max(Hin_arr);
            pert = 'perturbS';
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        elseif (FIG == 16)
            Hin_arr = [1.0:0.3:5.8, 5.8:-0.3:0.1, 0.1:0.3:1.0]*Hin_0;
            
            Hin_arr = Hin_arr(1:el_end);
            D0 = 1e-7*MORF*1;
            EndTime = 30000;
            NX = 256;
            NY = NX;
            el0 = 38; % First element of upward branch
            el1 = 17;  % Last element of upward branch
            LengthX = 50; % m
            DZ_TF = LengthX*slope*0.5; 
            DZ_MHWS = (MORF*Hin_arr) - DZ_TF;
            Hin_arr = -DZ_MHWS;
            xmin = min(Hin_arr); xmax = max(Hin_arr);
            pert = 'perturbS';
        elseif (FIG == 17)
            Hin_arr = [1.0:0.3:5.8, 5.8:-0.3:0.1, 0.1:0.3:1.0]*Hin_0;
            
            Hin_arr = Hin_arr(1:el_end);
            D0 = 1e-7*MORF*100;
            EndTime = 30000;
            NX = 256;
            NY = NX;
            el0 = 38; % First element of upward branch
            el1 = 17;  % Last element of upward branch
            LengthX = 50; % m
            DZ_TF = LengthX*slope*0.5; 
            DZ_MHWS = (MORF*Hin_arr) - DZ_TF;
            Hin_arr = -DZ_MHWS;
            xmin = min(Hin_arr); xmax = max(Hin_arr);
            pert = 'perturbS';
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        elseif (FIG >= 18) && (FIG <= 29) % NEED TO SPECIFY
            Hin_arr = [1.0:0.3:5.8, 5.8:-0.3:0.1, 0.1:0.3:1.0]*Hin_0;
            
            Hin_arr = Hin_arr(1:el_end);
            D0 = 0 %1e-7*MORF*100;
%             EndTime = 30000;
            NX = 256;
            NY = NX;
            el0 = 38; % First element of upward branch
            el1 = 17;  % Last element of upward branch
            LengthX = 50; % m
            DZ_TF = LengthX*slope*0.5; 
            DZ_MHWS = (MORF*Hin_arr) - DZ_TF;
            Hin_arr = -DZ_MHWS;
            xmin = min(Hin_arr); xmax = max(Hin_arr);
            pert = 'perturbS';
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        elseif (FIG >= 30) % NEED TO SPECIFY
            Hin_arr = [1.0:0.3:5.8, 5.8:-0.3:0.1, 0.1:0.3:1.0]*Hin_0;
            
            
            D0 = 0 %1e-7*MORF*100;
%             EndTime = 30000;
            NX = 256;
            NY = NX;
            el0 = 38; % First element of upward branch
            el1 = 17;  % Last element of upward branch
            LengthX = 50; % m
            DZ_TF = LengthX*slope*0.5; 
            DZ_MHWS = (MORF*Hin_arr) - DZ_TF;
            Hin_arr = -DZ_MHWS;
            xmin = min(Hin_arr); xmax = max(Hin_arr);
            Hin_arr = Hin_arr(1:el_end);
            pert = 'perturbS';
        end
        
        %%
        L = length(Hin_arr);
        Z = zeros(NX,NY,L);
        Z_mea = zeros(1,L);
        Z_min = zeros(1,L);
        Z_max = zeros(1,L);
        Z_var = zeros(1,L);
        Z_amp = zeros(1,L);
        Z_std = zeros(1,L);
        Z_qxy = zeros(1,L);
        Z_qxy_b0 = zeros(1,L);
        Z_qxy_v0 = zeros(1,L);
        Z_qxy_bv0 = zeros(1,L);
        Z_rr = zeros(1,L);
        
        %%
        for s_or_d = 1:2%1:2
            
            for SUBRUN = 0:L-1
                
                i = SUBRUN + 1;

                if s_or_d == 1
                    Ztemp = load(['s_FIG',num2str(FIG),'_SUBRUN',num2str(SUBRUN),'_SUBRUN2',num2str(SUBRUN2),'.mat']); Z(:,:,i) = Ztemp.s;
%                     Hin_order_i = [17:-1:1, 40:-1:37,    18:36];
                    
                    
%                     figure(444)
%                     subplot(5,9,i)
%                     Z0 = Z(:,:,i);
%                     imagesc(Z0);
%                     caxis([0, 0.3])
%                     set(gca,'xtick',[],'ytick',[])
%                     title(sprintf('%0.2f',Hin_arr(i)),'FontSize',15)
%                     axis image
                    
                elseif s_or_d == 2
                    Ztemp = load(['d_FIG',num2str(FIG),'_SUBRUN',num2str(SUBRUN),'_SUBRUN2',num2str(SUBRUN2),'.mat']); Z(:,:,i) = Ztemp.d;
                    
                    if ampvarmea == 4
                    figure(444)
                    subplot(5,9,i)
                    Z0 = Z(:,:,i);
                    imagesc(Z0);
                    caxis([0, 1500])
                    set(gca,'xtick',[],'ytick',[])
                    title(sprintf('%0.2f',Hin_arr(i)),'FontSize',15)
                    axis image
%                     if i > L-1
                        print(gcf,[titlestr,'_plot1.png'],'-dpng','-r500');
%                     end
                    end
                    
                    %%
                    Z1_order_i = [17:-1:1, 22:41,20:-1:18]
                    if i <= length(Z1_order_i)
                        if ampvarmea == 4
    %                         if i > 1
                            figure(555)
                            Z1 = Z(:,:,i);
                            subplot(2,21,Z1_order_i(i))
                            imagesc(Z1); axis image
                            title(sprintf('%0.2f',Hin_arr(i)),'FontSize',8)
                            set(gca,'xtick',[],'ytick',[])
                            caxis([0 1500])
    %                         if i > L-1
                            print(gcf,[titlestr,'_plot3.png'],'-dpng','-r500');
    %                         end
                                %                         saveas(gcf,[titlestr,'_plots2.png'])
    %                         end
                        end
                    end
                elseif s_or_d == 3
                    Ztemp = load(['d_FIG',num2str(FIG),'_SUBRUN',num2str(SUBRUN),'_SUBRUN2',num2str(SUBRUN2),'.mat']); Z(:,:,i) = Ztemp.d;
                    Ztemp = load(['h_FIG',num2str(FIG),'_SUBRUN',num2str(SUBRUN),'_SUBRUN2',num2str(SUBRUN2),'.mat']); Zh(:,:,i) = Ztemp.h;
                    Ztemp = load(['u_FIG',num2str(FIG),'_SUBRUN',num2str(SUBRUN),'_SUBRUN2',num2str(SUBRUN2),'.mat']); Zu(:,:,i) = Ztemp.u;
                    Ztemp = load(['v_FIG',num2str(FIG),'_SUBRUN',num2str(SUBRUN),'_SUBRUN2',num2str(SUBRUN2),'.mat']); Zv(:,:,i) = Ztemp.v;
                    Zd0 = Z(:,:,i);
                    Zh0 = Zh(:,:,i);
                    Zu0 = Zu(:,:,i);
                    Zv0 = Zv(:,:,i);
                    NS = sqrt(Zu0.*Zu0 + Zv0.*Zv0);
                    QXY = sum(Zh0(:).*NS(:))/sum(Zh0(:));
                    QXY = mean(Zv0(:));
                    RR = sum(Zh0(1,:).*Zv0(1,:))/sum(Zh0(:));
                    QXY_b0 = (sqrt(slope)/nb)*mean(Zh0(:)).^(1/6);
                    QXY_v0 = (sqrt(slope)/nv)*mean(Zh0(:)).^(1/6);
                    QXY_bv0 = (sqrt(slope)/(nb+(nv-nb)*(mean(Zd0(:))/kk)))*mean(Zh0(:)).^(1/6);
                end
                
                Z0 = Z(:,:,i);
                
                Z_mea(1,i) = mean(Z0(:));
                Z_min(1,i) = min(Z0(:));
                Z_max(1,i) = max(Z0(:));
                Z_std(1,i) = std(Z0(:));
                
%                 if ampvarmea == 5
%                 Z_qxy(1,i) = QXY;
%                 Z_rr(1,i) = RR;
%                 Z_qxy_b0(1,i) = QXY_b0;
%                 Z_qxy_v0(1,i) = QXY_v0;
%                 Z_qxy_bv0(1,i) = QXY_bv0;
%                 end
                
                [row col]=size(Z0);
                xlin=linspace(0,1,col);
                ylin=linspace(0,1,row);
                [X Y]=meshgrid(xlin,ylin);
                x=X(:);y=Y(:); z=Z0(:);
                xid=[ones(size(x)) y y.^2 y.^3 ]; %y.^4
%                 par=regress(z(:),xid);
%                 T=par(1)+par(2)*Y+par(3)*Y.^2+par(4)*Y.^3; %+par(5)*Y.^4;
%                 Zdet=Z0-T; Zdet = Zdet - min(Zdet(:));
%                 Z_var(1,i) = var(Zdet(:));
%                 Z_amp(1,i) = max(Zdet(:)) - min(Zdet(:));
                
                if plotting == 1
                    figure(i)
                    imagesc(Z0)
                    colorbar
                    axis image
                    if s_or_d == 1
                        caxis([0 1.1])
                        namestr = ['ART4_FIG',num2str(FIG),'_SUBRUN',num2str(SUBRUN),'_SUBRUN2',num2str(SUBRUN2),'_S'];
                    elseif s_or_d == 2
                        caxis([0 1500])
                        namestr = ['ART4_FIG',num2str(FIG),'_SUBRUN',num2str(SUBRUN),'_SUBRUN2',num2str(SUBRUN2),'_B'];
                    end
                    print('-dpng','-r300',namestr)
                end

            end % End SUBRUN-loop
        

        %% Sort data (up/downward; for shading purposes)
        if shade == 1
        elup = [el0:length(Hin_arr),1:el1];
        eldo = [el1+1:el0-1];
        Hin_arr_up = Hin_arr(elup);
        Hin_arr_do = Hin_arr(eldo);
        Z_mea_up = Z_mea(elup);
        Z_mea_do = Z_mea(eldo);
        Z_min_up = Z_min(elup);
        Z_min_do = Z_min(eldo);
        Z_max_up = Z_max(elup);
        Z_max_do = Z_max(eldo);
        Z_std_up = Z_std(elup);
        Z_std_do = Z_std(eldo);
        Z_qxy_up = Z_qxy(elup);
        Z_qxy_do = Z_qxy(eldo);
        Z_qxy_b0_up = Z_qxy_b0(elup);
        Z_qxy_b0_do = Z_qxy_b0(eldo);
        Z_qxy_v0_up = Z_qxy_v0(elup);
        Z_qxy_v0_do = Z_qxy_v0(eldo);
        Z_qxy_bv0_up = Z_qxy_bv0(elup);
        Z_qxy_bv0_do = Z_qxy_bv0(eldo);
        Z_rr_up = Z_rr(elup);
        Z_rr_do = Z_rr(eldo);
        end
        
        %%
        
        for ampvarmea = ampvarmea%1:3
            %%
            fs = 12;
            lw = 2;
            ms = 20;
            lw2 = 0.5;
            ms2 = 3;

            figure(101)
            xvec = linspace(0,3,2);
            yvec = zeros(1,length(xvec));
            x=xmin:0.1:xmax; y=0;

%             if (ampvarmea == 1) && (s_or_d) == 1
%                 subplot(2,3,1)
%                 if logplot == 1
%                     semilogx(Hin_arr,Z_amp,'.-k','LineWidth',lw,'MarkerSize',ms);
%                 elseif logplot == 0
%                     plot(Hin_arr,Z_amp,'.-k','LineWidth',lw,'MarkerSize',ms);
%                 end
%                 axis([xmin xmax 0 0.18])                                    % Samp
%                 ylabel('Detrended topographic amplitude [m]','FontSize',fs) % Samp

%             elseif (ampvarmea == 2) && (s_or_d) == 1
%                 subplot(2,3,2)
%                 if logplot == 1
%                     semilogx(Hin_arr,Z_var,'.-k','LineWidth',lw,'MarkerSize',ms);
%                 elseif logplot ==0
%                     plot(Hin_arr,Z_var,'.-k','LineWidth',lw,'MarkerSize',ms);
%                 end
%                 axis([xmin xmax 0 30e-4])                                   % Svar
%                 ylabel('Detrended topographic variance','FontSize',fs)      % Svar
                
                
            if (ampvarmea == 3) && (s_or_d) == 1
%                 subplot(2,3,3)
                subplot(2,1,1)
                if logplot == 1
                   
                    if shade == 0    
                        semilogx(Hin_arr,Z_mea,'.-b','LineWidth',lw,'MarkerSize',ms); hold on
                        semilogx(Hin_arr,Z_mea,'.-r','LineWidth',lw,'MarkerSize',ms);
                    elseif shade == 1
                        semilogx(Hin_arr_up,Z_mea_up,'.-b','LineWidth',lw,'MarkerSize',ms); hold on
                        semilogx(Hin_arr_do,Z_mea_do,'.-r','LineWidth',lw,'MarkerSize',ms);
                        
                        semilogx(Hin_arr_up, Z_min_up, '+-b', 'LineWidth', 1,'MarkerSize',ms2);%                     hold on;
                        semilogx(Hin_arr_up, Z_max_up, '+-b', 'LineWidth', 1,'MarkerSize',ms2);
                        x2 = [Hin_arr_up, fliplr(Hin_arr_up)];
                        inBetween = [Z_min_up, fliplr(Z_max_up)];
                        fill(x2, inBetween, 'b','FaceAlpha',0.15);

                        semilogx(Hin_arr_do, Z_min_do, 'x-r', 'LineWidth', 1,'MarkerSize',ms2);%                     hold on;
                        semilogx(Hin_arr_do, Z_max_do, 'x-r', 'LineWidth', 1,'MarkerSize',ms2);
                        x2 = [Hin_arr_do, fliplr(Hin_arr_do)];
                        inBetween = [Z_min_do, fliplr(Z_max_do)];
                        fill(x2, inBetween, 'r','FaceAlpha',0.15);
                    end
                    
                elseif logplot ==0
                    
                    if shade == 0
                        plot(Hin_arr,Z_mea,'.-b','LineWidth',lw,'MarkerSize',ms); hold on
                        S1_Hin = Hin_arr;
                        S1_mea = Z_mea;
                    elseif shade == 1
                        plot(Hin_arr_up,Z_mea_up,'.-b','LineWidth',lw,'MarkerSize',ms); hold on
                        plot(Hin_arr_do,Z_mea_do,'.-r','LineWidth',lw,'MarkerSize',ms);
                    
                        plot(Hin_arr_up, Z_min_up, '+-b', 'LineWidth', 1,'MarkerSize',ms2);%                     hold on;
                        plot(Hin_arr_up, Z_max_up, '+-b', 'LineWidth', 1,'MarkerSize',ms2);
                        x2 = [Hin_arr_up, fliplr(Hin_arr_up)];
                        inBetween = [Z_min_up, fliplr(Z_max_up)];
                        fill(x2, inBetween, 'b','FaceAlpha',0.15);

                        plot(Hin_arr_do, Z_min_do, 'x-r', 'LineWidth', 1,'MarkerSize',ms2);%                     hold on;
                        plot(Hin_arr_do, Z_max_do, 'x-r', 'LineWidth', 1,'MarkerSize',ms2);
                        x2 = [Hin_arr_do, fliplr(Hin_arr_do)];
                        inBetween = [Z_min_do, fliplr(Z_max_do)];
                        fill(x2, inBetween, 'r','FaceAlpha',0.15);
                    end
                   
                end
                axis([xmin xmax 0 0.2])                                   % Svar
                                axis([xmin xmax 0 1.1])                                   % Svar
                ylabel('Sediment elevation [m]','FontSize',fs)      % Svar
                legend({'\leftarrow (sinking)  ','\rightarrow (rising)  '},'FontSize',14,'Location','NorthWest')
                legend boxoff
    %                 xlabel('H_{in} [m/s]','FontSize',fs)
                set(gca,'FontSize',fs)
                title(['FIG',num2str(FIG),': D0 = 1e-7*MORF*',num2str(D0/MORF/(1e-7)),''])
                titlestr = ['ART4_FIG',num2str(FIG),'_SUBRUN2',num2str(SUBRUN2),'_bifurc_logplot=',num2str(logplot),'_abs'];

            elseif (ampvarmea == 3) && (s_or_d) == 2
%                 subplot(2,3,6)
                    subplot(2,1,2)
                if logplot == 1
                    
                     if shade == 0    
                        semilogx(Hin_arr,Z_mea,'.-b','LineWidth',lw,'MarkerSize',ms); hold on
                        semilogx(Hin_arr,Z_mea,'.-r','LineWidth',lw,'MarkerSize',ms);
                    elseif shade == 1
                        semilogx(Hin_arr_up,Z_mea_up,'.-b','LineWidth',lw,'MarkerSize',ms); hold on
                        semilogx(Hin_arr_do,Z_mea_do,'.-r','LineWidth',lw,'MarkerSize',ms);
                        
                        semilogx(Hin_arr_up, Z_min_up, '+-b', 'LineWidth', 1,'MarkerSize',ms2);%                     hold on;
                        semilogx(Hin_arr_up, Z_max_up, '+-b', 'LineWidth', 1,'MarkerSize',ms2);
                        x2 = [Hin_arr_up, fliplr(Hin_arr_up)];
                        inBetween = [Z_min_up, fliplr(Z_max_up)];
                        fill(x2, inBetween, 'b','FaceAlpha',0.15);

                        semilogx(Hin_arr_do, Z_min_do, 'x-r', 'LineWidth', 1,'MarkerSize',ms2);%                     hold on;
                        semilogx(Hin_arr_do, Z_max_do, 'x-r', 'LineWidth', 1,'MarkerSize',ms2);
                        x2 = [Hin_arr_do, fliplr(Hin_arr_do)];
                        inBetween = [Z_min_do, fliplr(Z_max_do)];
                        fill(x2, inBetween, 'r','FaceAlpha',0.15);
                    end
                    
                elseif logplot ==0
                    
                    if shade == 0
                        plot(Hin_arr,Z_mea,'.-b','LineWidth',lw,'MarkerSize',ms);
                    elseif shade == 1
                        kk = 1500/100;
                        plot(Hin_arr_up,Z_mea_up/kk,'.-b','LineWidth',lw,'MarkerSize',ms); hold on
                        plot(Hin_arr_do,Z_mea_do/kk,'.-r','LineWidth',lw,'MarkerSize',ms);
                    
                        plot(Hin_arr_up, Z_min_up/kk, '+-b', 'LineWidth', lw2,'MarkerSize',ms2);%                     hold on;
                        plot(Hin_arr_up, Z_max_up/kk, '+-b', 'LineWidth', lw2,'MarkerSize',ms2);
                        x2 = [Hin_arr_up, fliplr(Hin_arr_up)];
                        inBetween = [Z_min_up/kk, fliplr(Z_max_up/kk)];
                        fill(x2, inBetween, 'b','FaceAlpha',0.15);

                        plot(Hin_arr_do, Z_min_do/kk, 'x-r', 'LineWidth', lw2,'MarkerSize',ms2);%                     hold on;
                        plot(Hin_arr_do, Z_max_do/kk, 'x-r', 'LineWidth', lw2,'MarkerSize',ms2);
                        x2 = [Hin_arr_do, fliplr(Hin_arr_do)];
                        inBetween = [Z_min_do/kk, fliplr(Z_max_do/kk)];
                        fill(x2, inBetween, 'r','FaceAlpha',0.15);
                    end

                end
                
                axis([xmin xmax 0 100])                                    % Bvar
                ylabel('Algal density [% of carr. cap.]','FontSize',fs)          % Bvar
%                 legend({'{\leftarrow} (sinking)  ','\rightarrow (rising)  '},'FontSize',14,'Location','NorthWest')
%                 legend boxoff
                xlabel('Marsh elevation relative to MHWS [m]','FontSize',fs)
                set(gca,'FontSize',fs)
%                 title(['FIG',num2str(FIG),': D0 = 1e-7*MORF*',num2str(D0/MORF/(1e-7)),''])
            titlestr = ['ART4_FIG',num2str(FIG),'_SUBRUN2',num2str(SUBRUN2),'_bifurc_logplot=',num2str(logplot),'_abs'];


















            elseif (ampvarmea == 4) && (s_or_d) == 1
                subplot(2,1,1)
                
                plot(Hin_arr_up,Z_std_up,'.-b','LineWidth',lw,'MarkerSize',ms); hold on
                plot(Hin_arr_do,Z_std_do,'.-r','LineWidth',lw,'MarkerSize',ms);

                axis([xmin xmax 0 0.10])                                   % Svar
                ylabel('Sediment elevation std [m]','FontSize',fs)      % Svar
                legend({'\leftarrow (sinking)  ','\rightarrow (rising)  '},'FontSize',14,'Location','NorthWest')
                legend boxoff
    %                 xlabel('H_{in} [m/s]','FontSize',fs)
                set(gca,'FontSize',fs)
                title(['FIG',num2str(FIG),': D0 = 1e-7*MORF*',num2str(D0/MORF/(1e-7)),''])
                titlestr = ['ART4_FIG',num2str(FIG),'_SUBRUN2',num2str(SUBRUN2),'_bifurc_logplot=',num2str(logplot),'_std'];

            elseif (ampvarmea == 4) && (s_or_d) == 2
                subplot(2,1,2)
                    
                plot(Hin_arr_up,Z_std_up,'.-b','LineWidth',lw,'MarkerSize',ms); hold on
                plot(Hin_arr_do,Z_std_do,'.-r','LineWidth',lw,'MarkerSize',ms);
                
                axis([xmin xmax 0 600])                                    % Bvar
                ylabel('Vegetation density std [m^{-2}]','FontSize',fs)          % Bvar
                legend({'{\leftarrow} (sinking)  ','\rightarrow (rising)  '},'FontSize',14,'Location','NorthWest')
                legend boxoff
                xlabel('Marsh elevation relative to MHWS [m]','FontSize',fs)
                set(gca,'FontSize',fs)
%                 title(['FIG',num2str(FIG),': D0 = 1e-7*MORF*',num2str(D0/MORF/(1e-7)),''])
                titlestr = ['ART4_FIG',num2str(FIG),'_SUBRUN2',num2str(SUBRUN2),'_bifurc_logplot=',num2str(logplot),'_std'];
                
            elseif (ampvarmea == 5) && (s_or_d == 3)
                subplot(2,1,1)
                plot(Hin_arr_up,Z_rr_up,'.-b','LineWidth',lw,'MarkerSize',ms); hold on
                plot(Hin_arr_do,Z_rr_do,'.-r','LineWidth',lw,'MarkerSize',ms);
                axis([xmin xmax 0 2.5e-3])                                    % Bvar
                ylabel('Refreshment rate [1/s]','FontSize',fs)          % Bvar
                legend({'{\leftarrow} (sinking)  ','\rightarrow (rising)  '},'FontSize',14,'Location','NorthEast')
                legend boxoff
%                 xlabel('Marsh elevation relative to MHWS [m]','FontSize',fs)
                set(gca,'FontSize',fs)
%                 title(['FIG',num2str(FIG),': D0 = 1e-7*MORF*',num2str(D0/MORF/(1e-7)),''])
                            pbaspect([2 1 1])

                subplot(2,1,2)
                plot(Hin_arr_up,Z_qxy_up,'.-b','LineWidth',lw,'MarkerSize',ms); hold on
                plot(Hin_arr_do,Z_qxy_do,'.-r','LineWidth',lw,'MarkerSize',ms);
%                 plot(Hin_arr_up,Z_qxy_b0_up,'-k','LineWidth',lw,'MarkerSize',ms);
%                 plot(Hin_arr_do,Z_qxy_b0_do,'-k','LineWidth',lw,'MarkerSize',ms);
%                 plot(Hin_arr_up,Z_qxy_v0_up,'-k','LineWidth',lw,'MarkerSize',ms);
%                 plot(Hin_arr_do,Z_qxy_v0_do,'-k','LineWidth',lw,'MarkerSize',ms);
                plot(Hin_arr_up,Z_qxy_bv0_up,'-k','LineWidth',lw,'MarkerSize',ms);
                plot(Hin_arr_do,Z_qxy_bv0_do,'-k','LineWidth',lw,'MarkerSize',ms);
%                 axis([xmin xmax 0 0.8])                                    % Bvar
                ylabel('Drainage flow speed [m/s]','FontSize',fs)          % Bvar
                legend({'{\leftarrow} (sinking)  ','\rightarrow (rising)  '},'FontSize',14,'Location','NorthEast')
                legend boxoff
                xlabel('Marsh elevation relative to MHWS [m]','FontSize',fs)
                set(gca,'FontSize',fs)
                title(['FIG',num2str(FIG),': D0 = 1e-7*MORF*',num2str(D0/MORF/(1e-7)),''])
                titlestr = ['ART4_FIG',num2str(FIG),'_SUBRUN2',num2str(SUBRUN2),'_bifurc_logplot=',num2str(logplot),'_qxy'];
                            pbaspect([2 1 1])

            end    
            pbaspect([2 1 1])

        end
        
    end
        
print('-dpng','-r300',titlestr)
% print('-dpdf','-r300',titlestr)

%%
if ampvarmea == 3
    figure(999)
    nrow=5;
    ncol=5;
    subplot(nrow,ncol,SUBRUN2)               
    plot(Hin_arr_up,Z_mea_up,'.-b','LineWidth',lw,'MarkerSize',ms); hold on
                        plot(Hin_arr_do,Z_mea_do,'.-r','LineWidth',lw,'MarkerSize',ms);
                    
                        plot(Hin_arr_up, Z_min_up, '+-b', 'LineWidth', lw2,'MarkerSize',ms2);%                     hold on;
                        plot(Hin_arr_up, Z_max_up, '+-b', 'LineWidth', lw2,'MarkerSize',ms2);
                        x2 = [Hin_arr_up, fliplr(Hin_arr_up)];
                        inBetween = [Z_min_up, fliplr(Z_max_up)];
                        fill(x2, inBetween, 'b','FaceAlpha',0.15);

                        plot(Hin_arr_do, Z_min_do, 'x-r', 'LineWidth', lw2,'MarkerSize',ms2);%                     hold on;
                        plot(Hin_arr_do, Z_max_do, 'x-r', 'LineWidth', lw2,'MarkerSize',ms2);
                        x2 = [Hin_arr_do, fliplr(Hin_arr_do)];
                        inBetween = [Z_min_do, fliplr(Z_max_do)];
                        fill(x2, inBetween, 'r','FaceAlpha',0.15);
                        
    axis([xmin xmax 0 1500])                                    % Bvar
    % ylabel('Vegetation density std [m^{-2}]','FontSize',fs)          % Bvar
    % legend({'{\leftarrow} (sinking)  ','\rightarrow (rising)  '},'FontSize',14,'Location','NorthWest')
    % legend boxoff
    % xlabel('Marsh elevation relative to MHWS [m]','FontSize',fs)
    title(SUBRUN2,'FontSize',10)
    set(gca,'FontSize',fs,'XTick',[],'YTick',[])
    %                 title(['FIG',num2str(FIG),': D0 = 1e-7*MORF*',num2str(D0/MORF/(1e-7)),''])
    titlestr = ['ART4_FIG',num2str(FIG),'_SUBRUN2',num2str(SUBRUN2),'_bifurc_logplot=',num2str(logplot),'_abs_overview'];
    saveas(gcf,[titlestr,'.png'])
elseif ampvarmea == 4
    figure(999)
    nrow=6;
    ncol=5;
    subplot(nrow,ncol,SUBRUN2)               
    plot(Hin_arr_up,Z_std_up,'.-b','LineWidth',lw,'MarkerSize',ms); hold on
    plot(Hin_arr_do,Z_std_do,'.-r','LineWidth',lw,'MarkerSize',ms);

    axis([xmin xmax 0 800])                                    % Bvar
    % ylabel('Vegetation density std [m^{-2}]','FontSize',fs)          % Bvar
    % legend({'{\leftarrow} (sinking)  ','\rightarrow (rising)  '},'FontSize',14,'Location','NorthWest')
    % legend boxoff
    % xlabel('Marsh elevation relative to MHWS [m]','FontSize',fs)
    title(SUBRUN2,'FontSize',10)
    set(gca,'FontSize',fs,'XTick',[],'YTick',[])
    %                 title(['FIG',num2str(FIG),': D0 = 1e-7*MORF*',num2str(D0/MORF/(1e-7)),''])
    titlestr = ['ART4_FIG',num2str(FIG),'_SUBRUN2',num2str(SUBRUN2),'_bifurc_logplot=',num2str(logplot),'_std_overview'];
    saveas(gcf,[titlestr,'.png'])
end
end


% %%
% % i = 1;
% MORF = 44712;
% Hin_arr = [1.0:0.3:5.8, 5.8:-0.3:0.1, 0.1:0.3:1.0]*1.4e-5;
% 
% Df = 1.4e-5/0.0225;
% nb = 0.016;
% nv = 0.2;
% kk = 1500;
% S_in = MORF*5e-9;
% Hc = 0.001;
% Qs = 6e-4;
% Es = MORF*3.5e-4;
% E0 = Es;
% pE = 0.9;
% g = 9.81;
% Qq = 0.16;
% Eb = MORF*2e-5;
% Ec = Eb;
% rr = MORF*3.2e-8;
% 
% h0_arr = Hin_arr/Df;
% u0_arr = (sqrt(slope)/nb)*h0_arr.^(1/6);
% S0_arr = (S_in/(Es*g*slope))*(h0_arr-Hc)./(Qs+h0_arr-Hc);
% 
% Hin_arr = -DZ_MHWS;
% % Hin_arr = 
% 
% figure
% subplot(2,2,1)
% plot(Hin_arr,u0_arr,'b*-')
% title('u0')
% 
% subplot(2,2,2)
% plot(Hin_arr,h0_arr,'b*-')
% title('h0')
% 
% subplot(2,2,3)
% plot(Hin_arr,S0_arr,'b*-')
% title('S0')
% 
% % subplot(2,2,1)
% % plot(Hin_arr,u0_arr,'b*-')
% % title('u0')

%%

% for i = 1:length(Hin_arr)
%     
%     Hin = Hin_arr(i);
%     
%     h0 = Hin/Df;
% %     syms UU_HOM0 SS_HOM0 BB_HOM0
% %     HomogenSolve0 = vpasolve([...    % Solve system of 2 equations numerically
% %         (UU_HOM0 == (sqrt(slope)./(nb+(nv-nb).*0./kk)).*h0.^(1/6)), ...
% %         (SS_HOM0 == ((S_in.*(h0-Hc)./(Qs+h0-Hc)).*h0.^(1/3)./(nb+(nv-nb).*0./kk).^2)./(Es.*(1-pE.*0./kk).*g.*UU_HOM0.^2))], ...
% %         [UU_HOM0, SS_HOM0],[1e-4,0.1])    ;   % OK!
% %     
% %     u00(i) = HomogenSolve0.UU_HOM0;
% %     v00(i) = 0;
% %     h00(i) = h0;
% %     S00(i) = HomogenSolve0.SS_HOM0;
% %     B00(i) = 0; %HomogenSolve0.BB_HOM0;
%     
% %     syms UU_HOM1 SS_HOM1 BB_HOM1 
% %     HomogenSolve1 = vpasolve([...    % Solve system of 2 equations numerically
% % %        N_HOM = (nb+(nv-nb)*BB_HOM1/kk);
% %         UU_HOM1 == (sqrt(slope)./(nb+(nv-nb).*BB_HOM1./kk)).*h0.^(1/6), ...
% %         SS_HOM1 == ((S_in.*(h0-Hc)./(Qs+h0-Hc)).*h0.^(1/3)./(nb+(nv-nb).*BB_HOM1./kk).^2)./(Es.*(1-pE.*BB_HOM1./kk).*g.*UU_HOM1.^2), ...
% %         BB_HOM1 == kk.*(1-(Qq+h0-Hc).*(Eb.*g.*UU_HOM1.^2)./(rr.*Qq.*h0.^(1/3)./(nb+(nv-nb).*BB_HOM1./kk).^2)) ], ...
% %         [UU_HOM1, SS_HOM1, BB_HOM1],[0.5,0.1,1000])  ;     % OK!
% %    
% %     u01(i) = HomogenSolve1.UU_HOM1;
% %     v01(i) = 0;
% %     h01(i) = h0;
% %     S01(i) = HomogenSolve1.SS_HOM1;
% %     B01(i) = HomogenSolve1.BB_HOM1;
% 
% %     BB_HOM1 = 0;
% %     syms UU_HOM1 SS_HOM1
% %     HomogenSolve1 = vpasolve([...    % Solve system of 2 equations numerically
% %         UU_HOM1 == (sqrt(slope)./(nb+(nv-nb).*BB_HOM1./kk)).*h0.^(1/6), ...
% %         0 == S_in*(h0-Hc)/(Qs+h0-Hc) - E0*(1.0-pE*BB_HOM1/kk)*SS_HOM1*UU_HOM1*UU_HOM1*(g/((1/(nb+(nv-nb)*BB_HOM1/kk))*h0^(1/6))^2)], ...
% %         [UU_HOM1, SS_HOM1])  ;     % OK!
% 
%     syms UU_HOM1 SS_HOM1 BB_HOM1 
%     HomogenSolve1 = vpasolve([...    % Solve system of 2 equations numerically
%         UU_HOM1 == (sqrt(slope)./(nb+(nv-nb).*BB_HOM1./kk)).*h0.^(1/6), ...
%         0 == S_in*(h0-Hc)/(Qs+h0-Hc) - E0*(1.0-pE*BB_HOM1/kk)*SS_HOM1*UU_HOM1*UU_HOM1*(g/((1/(nb+(nv-nb)*BB_HOM1/kk))*h0^(1/6))^2), ...
%         0 == rr*BB_HOM1*(1 - BB_HOM1/kk)*Qq/(Qq+h0-Hc) - Ec*BB_HOM1*UU_HOM1*UU_HOM1*(g/((1/(nb+(nv-nb)*BB_HOM1/kk))*h0^(1/6))^2)], ...
%         [UU_HOM1, SS_HOM1, BB_HOM1],[1,0.1,3])  ;     % OK!
%    
%     u00(i) = HomogenSolve1.UU_HOM1(1);
%     v00(i) = 0;
%     h00(i) = h0;
%     S00(i) = HomogenSolve1.SS_HOM1(1);
%     B00(i) = HomogenSolve1.BB_HOM1(1);
%     
%     u01(i) = HomogenSolve1.UU_HOM1(2);
%     v01(i) = 0;
%     h01(i) = h0;
%     S01(i) = HomogenSolve1.SS_HOM1(2);
%     B01(i) = HomogenSolve1.BB_HOM1(2);
% end
% %%
% close(figure(999));figure(999)
% Hin_arr = -DZ_MHWS;
% subplot(2,2,1)
% plot(Hin_arr,u00,'-*k'); hold on
% % plot(Hin_arr,u01,'-*g');
% title('u')
% 
% subplot(2,2,2)
% plot(Hin_arr,h00,'-*k'); hold on
% % plot(Hin_arr,h01,'-*g');
% title('h')
% 
% subplot(2,2,3)
% plot(Hin_arr,S00,'-*k'); hold on
% % plot(Hin_arr,S01,'-*g');
% title('S')
% 
% subplot(2,2,4)
% plot(Hin_arr,B00,'-*k'); hold on
% % plot(Hin_arr,B01,'-*g');
% title('B')
% saveas(gcf,'HomogeneousSolution.png')

%%
% -DZ_MHWS
% close(figure(44));figure(44)
% fs = 18
% col = [0.5 0.5 0.5];
% x1 = [-0.041,-0.7825,-1.9093,-2.2848,-2.4726];
% x2 = [-0.0313,-0.5947,-1.3459,-1.3459,-1.3459];
% y = [1000,100,10,1,1/10];
% 
% semilogy([x1(1),x2(1)],[1000,1000],'-','Color',col,'LineWidth',10); hold on
% semilogy([x1(2),x2(2)],[100,100],'-','Color',col,'LineWidth',10)
% semilogy([x1(3),x2(3)],[10,10],'-','Color',col,'LineWidth',10)
% semilogy([x1(4),x2(4)],[1,1],'-','Color',col,'LineWidth',10)
% semilogy([x1(5),x2(5)],[1/10,1/10],'-','Color',col,'LineWidth',10);
% semilogy(x1,y,'g-')
% semilogy(x2,y,'g-')
% axis([xmin xmax 1/20 2000])
% xlabel('Marsh elevation relative to MHWS [m]','FontSize',fs)
% ylabel('Topographic diffusivity, D0 [m^2/s]','FontSize',fs)
% ll = legend({['Bistability' char(10) 'range']},'Location','SouthEast')
% set(ll,'FontSize',14)
% % legend boxoff
% set(gca,'FontSize',fs)
% pbaspect([2 1 1])
% % saveas(gcf,'RangeOfBistability.png')
% print('-dpng','-r300','RangeOfBistability_pD=0')
end
