% Compare maximum cross-sectional flow above baseflow

clear

max_dist = 27000;
min_t = 7200; % plots start at 2 hour
max_t = 96 * 3600; % hours * seconds in an hour - gives seconds
peak_index = [180; 300; 180; 300; 180; 300];

% Load cross-section ids
id_floc = '/Users/Colin/Documents/MATLAB/DFlow/influences_chap2/';
load(strcat(id_floc,'xs_xs_ids.mat'));
load(strcat(id_floc,'xs_rcl_ids.mat'));

%% FP n = 0.040
% Collect baseflow maximum for FP n = 0.040

base30cms040_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3095_30cmsfp040/'...
    'abq10a12bf30cmsfp040rst01_0000_20160104_040000_his.nc'];
base150cms040_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3096_150cmsfp040/'...
    'abq10a12bf150cmsfp040rst01_0000_20160103_000000_his.nc'];

% Extract flow at all cross-section times
base30cms040_all_Q = ncread(base30cms040_file,'cross_section_discharge');
base150cms040_all_Q = ncread(base150cms040_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
base30cms040_Q = base30cms040_all_Q(:,end);
base150cms040_Q = base150cms040_all_Q(:,end);

clear base30cms040_all_Q base150cms040_all_Q

%% Collect flow under steady peak flow for FP n = 0.040

peak180cms040_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3107_f180cmsfp040/'...
    'abq10a12bf180cmsfp040rst02_0000_20160104_130000_his.nc'];
peak300cms040_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3108_f300cmsfp040/'...
    'abq10a12bf300cmsfp040rst02_0000_20160103_230000_his.nc'];

% Extract flow at all cross-section times
peak180cms040_all_Q = ncread(peak180cms040_file,'cross_section_discharge');
peak300cms040_all_Q = ncread(peak300cms040_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
peak180cms040_Q = peak180cms040_all_Q(:,end);
peak300cms040_Q = peak300cms040_all_Q(:,end);

clear peak180cms040_all_Q peak300cms040_all_Q

%% Collect peak flow information from FP n = 0.040 flood runs

fpath_root = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/'];

% Define folder with *his.nc files
fpath = ['run3101_s180fp040/';'run3102_s300fp040/';...
    'run3103_m180fp040/';'run3104_m300fp040/';...
    'run3105_l180fp040/';'run3106_l300fp040/']; % file path

folder_num = numel(fpath(:,1));
xs_fp040_cell = cell(folder_num,7);

for fi = 1:folder_num
    
    xs_fp040_cell(fi,1) = {fpath(fi,:)};
    upstream_peak = peak_index(fi);

    if mod(fi,2) == 0
        xs_peak_Q = peak300cms040_Q;
        xs_base_Q = base150cms040_Q;
    elseif mod(fi,2) == 1
        xs_peak_Q = peak180cms040_Q;
        xs_base_Q = base30cms040_Q;
    end
    
    fpath_fi = strcat(fpath_root,fpath(fi,:));
    
    fnames = struct2cell(dir(strcat(fpath_fi,'*his.nc')))'; 
    [fnum,~] = size(fnames); % number of files

% Extract cross-sectional information

    % Returns cell array structured as the following
    % xs_cell
    % Column 1 - name of model run (assigned earlier)
    % Column 2 - Model time
    % Column 3 - XS x and y coordinates
    % Column 4 - All XS measured flow at all times
    % Column 5 - XS #, Time of max discharge, Maximum discharge
    % Column 6 - XS distance
    % Column 7 - Fraction of upstream peak discharge
    % Column 8 - Flow above baseflow
    
    [xs_fp040_cell(fi,2:10), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his_chap2(fnum, fnames, fpath_fi, max_dist, min_t,...
        max_t, xs_peak_Q, xs_base_Q, upstream_peak);
    
end

%% FP n = 0.060
% Collect baseflow maximum for FP n = 0.060

base30cms060_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3047_30cmsbank060/'...
    'abq10a12bf30cmsfp060rst01_0000_20160103_220000_his.nc'];
base150cms060_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3048_150cmsbank060/'...
    'abq10a12bf150cmsfp060rst01_0000_20160103_190000_his.nc'];

% Extract flow at all cross-section times
base30cms060_all_Q = ncread(base30cms060_file,'cross_section_discharge');
base150cms060_all_Q = ncread(base150cms060_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
base30cms060_Q = base30cms060_all_Q(:,end);
base150cms060_Q = base150cms060_all_Q(:,end);

clear base30cms060_all_Q base150cms060_all_Q

%% Collect flow under steady peak flow for FP n = 0.060

peak180cms060_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3049_180cmsbank060/'...
    'abq10a12b180cmsrst02_0000_20160104_140000_his.nc'];
peak300cms060_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3050_300cmsbank060/'...
    'abq10a12bf300cmsfp060rst02_0000_20160103_200000_his.nc'];

% Extract flow at all cross-section times
peak180cms060_all_Q = ncread(peak180cms060_file,'cross_section_discharge');
peak300cms060_all_Q = ncread(peak300cms060_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
peak180cms060_Q = peak180cms060_all_Q(:,end);
peak300cms060_Q = peak300cms060_all_Q(:,end);

clear peak180cms060_all_Q peak300cms060_all_Q

%% Collect peak flow information from FP n = 0.060 flood runs

fpath_root = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/'];

% Define folder with *his.nc files
fpath = ['run3053_s180fp060/';'run3054_s300fp060/';...
    'run3055_m180fp060/';'run3056_m300fp060/';...
    'run3057_l180fp060/';'run3058_l300fp060/']; % file path

folder_num = numel(fpath(:,1));
xs_fp060_cell = cell(folder_num,7);

for fi = 1:folder_num
    
    xs_fp060_cell(fi,1) = {fpath(fi,:)};
    upstream_peak = peak_index(fi);

    if mod(fi,2) == 0
        xs_peak_Q = peak300cms060_Q;
        xs_base_Q = base150cms060_Q;
    elseif mod(fi,2) == 1
        xs_peak_Q = peak180cms060_Q;
        xs_base_Q = base30cms060_Q;
    end
    
    fpath_fi = strcat(fpath_root,fpath(fi,:));
    
    fnames = struct2cell(dir(strcat(fpath_fi,'*his.nc')))'; 
    [fnum,~] = size(fnames); % number of files

% Extract cross-sectional information
    % Returns cell array structured as the following
    % xs_cell
    % Column 1 - name of model run (assigned earlier)
    % Column 2 - Model time
    % Column 3 - XS x and y coordinates
    % Column 4 - All XS measured flow at all times
    % Column 5 - XS #, Time of max discharge, Maximum discharge
    % Column 6 - XS distance
    % Column 7 - Fraction of upstream peak discharge
    % Column 8 - Flow above baseflow
    
    [xs_fp060_cell(fi,2:10), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his_chap2(fnum, fnames, fpath_fi, max_dist, min_t,...
        max_t, xs_peak_Q, xs_base_Q, upstream_peak);
    
end

%% FP n = 0.080
% Collect baseflow maximum for FP n = 0.080

base30cms080_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3097_30cmsfp080/'...
    'abq10a12bf30cmsfp080rst01_0000_20160104_070000_his.nc'];
base150cms080_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3098_150cmsfp080/'...
    'abq10a12bf150cmsfp080rst01_0000_20160103_000000_his.nc'];

% Extract flow at all cross-section times
base30cms080_all_Q = ncread(base30cms080_file,'cross_section_discharge');
base150cms080_all_Q = ncread(base150cms080_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
base30cms080_Q = base30cms080_all_Q(:,end);
base150cms080_Q = base150cms080_all_Q(:,end);

clear base30cms080_all_Q base150cms080_all_Q

%% Collect flow under steady peak flow for FP n = 0.080

peak180cms080_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3109_f180cmsfp080/'...
    'abq10a12bf180cmsfp080rst02_0000_20160104_120000_his.nc'];
peak300cms080_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3110_f300cmsfp080/'...
    'abq10a12bf300cmsfp080rst02_0000_20160103_200000_his.nc'];

% Extract flow at all cross-section times
peak180cms080_all_Q = ncread(peak180cms080_file,'cross_section_discharge');
peak300cms080_all_Q = ncread(peak300cms080_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
peak180cms080_Q = peak180cms080_all_Q(:,end);
peak300cms080_Q = peak300cms080_all_Q(:,end);

clear peak180cms080_all_Q peak300cms080_all_Q

%% Collect peak flow information from FP n = 0.080 flood runs

fpath_root = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/'];

% Define folder with *his.nc files
fpath = ['run3113_s180fp080/';'run3114_s300fp080/';...
    'run3115_m180fp080/';'run3116_m300fp080/';...
    'run3117_l180fp080/';'run3118_l300fp080/']; % file path

folder_num = numel(fpath(:,1));
xs_fp080_cell = cell(folder_num,7);

for fi = 1:folder_num
    
    xs_fp080_cell(fi,1) = {fpath(fi,:)};
    upstream_peak = peak_index(fi);

    if mod(fi,2) == 0
        xs_peak_Q = peak300cms080_Q;
        xs_base_Q = base150cms080_Q;
    elseif mod(fi,2) == 1
        xs_peak_Q = peak180cms080_Q;
        xs_base_Q = base30cms080_Q;
    end
    
    fpath_fi = strcat(fpath_root,fpath(fi,:));
    
    fnames = struct2cell(dir(strcat(fpath_fi,'*his.nc')))'; 
    [fnum,~] = size(fnames); % number of files

% Extract cross-sectional information
    % Returns cell array structured as the following
    % xs_cell
    % Column 1 - name of model run (assigned earlier)
    % Column 2 - Model time
    % Column 3 - XS x and y coordinates
    % Column 4 - All XS measured flow at all times
    % Column 5 - XS #, Time of max discharge, Maximum discharge
    % Column 6 - XS distance
    % Column 7 - Fraction of upstream peak discharge
    % Column 8 - Flow above baseflow
    
    [xs_fp080_cell(fi,2:10), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his_chap2(fnum, fnames, fpath_fi, max_dist, min_t,...
        max_t, xs_peak_Q, xs_base_Q, upstream_peak);
    
end

%% FP n = 0.100
% Collect baseflow maximum for FP n = 0.100

base30cms100_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3099_30cmsfp100/'...
    'abq10a12bf30cmsfp100rst01_0000_20160104_100000_his.nc'];
base150cms100_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3100_150cmsfp100/'...
    'abq10a12bf150cmsfp100rst02_0000_20160104_210000_his.nc'];

% Extract flow at all cross-section times
base30cms100_all_Q = ncread(base30cms100_file,'cross_section_discharge');
base150cms100_all_Q = ncread(base150cms100_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
base30cms100_Q = base30cms100_all_Q(:,end);
base150cms100_Q = base150cms100_all_Q(:,end);

clear base30cms100_all_Q base150cms100_all_Q

%% Collect flow under steady peak flow for FP n = 0.100

peak180cms100_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3111_f180cmsfp100/'...
    'abq10a12bf180cmsfp100rst02_0000_20160104_110000_his.nc'];
peak300cms100_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3112_f300cmsfp100/'...
    'abq10a12bf300cmsfp100rst03_0000_20160104_150000_his.nc'];

% Extract flow at all cross-section times
peak180cms100_all_Q = ncread(peak180cms100_file,'cross_section_discharge');
peak300cms100_all_Q = ncread(peak300cms100_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
peak180cms100_Q = peak180cms100_all_Q(:,end);
peak300cms100_Q = peak300cms100_all_Q(:,end);

clear peak180cms080_all_Q peak300cms080_all_Q

%% Collect peak flow information from FP n = 0.100 flood runs

fpath_root = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/'];

% Define folder with *his.nc files
fpath = ['run3119_s180fp100/';'run3120_s300fp100/';...
    'run3121_m180fp100/';'run3122_m300fp100/';...
    'run3123_l180fp100/';'run3124_l300fp100/']; % file path

folder_num = numel(fpath(:,1));
xs_fp100_cell = cell(folder_num,7);

for fi = 1:folder_num
    
    xs_fp100_cell(fi,1) = {fpath(fi,:)};
    upstream_peak = peak_index(fi);

    if mod(fi,2) == 0
        xs_peak_Q = peak300cms100_Q;
        xs_base_Q = base150cms100_Q;
    elseif mod(fi,2) == 1
        xs_peak_Q = peak180cms100_Q;
        xs_base_Q = base30cms100_Q;
    end
    
    fpath_fi = strcat(fpath_root,fpath(fi,:));
    
    fnames = struct2cell(dir(strcat(fpath_fi,'*his.nc')))'; 
    [fnum,~] = size(fnames); % number of files

% Extract cross-sectional information
    % Returns cell array structured as the following
    % xs_cell
    % Column 1 - name of model run (assigned earlier)
    % Column 2 - Model time
    % Column 3 - XS x and y coordinates
    % Column 4 - All XS measured flow at all times
    % Column 5 - XS #, Time of max discharge, Maximum discharge
    % Column 6 - XS distance
    % Column 7 - Fraction of upstream peak discharge
    % Column 8 - Flow above baseflow
    
    [xs_fp100_cell(fi,2:10), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his_chap2(fnum, fnames, fpath_fi, max_dist, min_t,...
        max_t, xs_peak_Q, xs_base_Q, upstream_peak);
    
end

%% Plot total floodplain flow
figure(1)
set(gcf,'color','white');
xs_dist = (250:250:max_dist)/1000;
for i = 1:6
    
    all_time = xs_fp040_cell{:,2};

    above_base_hold = xs_fp040_cell{i,4}';
    rfp_flow = above_base_hold(xs_rcl_ids(1:324,1)==1,:);
    lfp_flow = above_base_hold(xs_rcl_ids(1:324,1)==3,:);
    fp040_flow = rfp_flow + lfp_flow;
    ch_flow040 = above_base_hold(xs_rcl_ids(1:324,1)==2,:);
    clear above_base_hold rfp_flow lfp_flow
    peak_idtQ_040 = xs_fp040_cell{i,8};
    peak_t_040 = peak_idtQ_040(:,2);
    for j = 1:max_dist/250
        [peak_ti_040(j,1),~] = find(peak_t_040(j,1)==all_time(:,1));
        peak_fp040(j,1) = fp040_flow(j,peak_ti_040(j,1));
        peak_ch040(j,1) = ch_flow040(j,peak_ti_040(j,1));
        fp2ch_frac040(j,1) = peak_fp040(j,1) ./ peak_ch040(j,1);
        fp2tot_frac040(j,1) = peak_fp040(j,1) ./ ...
            (peak_ch040(j,1) + peak_fp040(j,1));
    end
    
    above_base_hold = xs_fp060_cell{i,4}';
    rfp_flow = above_base_hold(xs_rcl_ids(1:324,1)==1,:);
    lfp_flow = above_base_hold(xs_rcl_ids(1:324,1)==3,:);
    fp060_flow = rfp_flow + lfp_flow;
    ch_flow060 = above_base_hold(xs_rcl_ids(1:324,1)==2,:);
    clear above_base_hold rfp_flow lfp_flow
    peak_idtQ_060 = xs_fp060_cell{i,8};
    peak_t_060 = peak_idtQ_060(:,2);
    for j = 1:max_dist/250
        [peak_ti_060(j,1),~] = find(peak_t_060(j,1)==all_time(:,1));
        peak_fp060(j,1) = fp060_flow(j,peak_ti_060(j,1));
        peak_ch060(j,1) = ch_flow060(j,peak_ti_060(j,1));
        fp2ch_frac060(j,1) = peak_fp060(j,1) ./ peak_ch060(j,1);
        fp2tot_frac060(j,1) = peak_fp060(j,1) ./ ...
            (peak_ch060(j,1) + peak_fp060(j,1));
    end
    
    above_base_hold = xs_fp080_cell{i,4}';
    rfp_flow = above_base_hold(xs_rcl_ids(1:324,1)==1,:);
    lfp_flow = above_base_hold(xs_rcl_ids(1:324,1)==3,:);
    fp080_flow = rfp_flow + lfp_flow;
    ch_flow080 = above_base_hold(xs_rcl_ids(1:324,1)==2,:);
    clear above_base_hold rfp_flow lfp_flow
    peak_idtQ_080 = xs_fp080_cell{i,8};
    peak_t_080 = peak_idtQ_080(:,2);
    for j = 1:max_dist/250
        [peak_ti_080(j,1),~] = find(peak_t_080(j,1)==all_time(:,1));
        peak_fp080(j,1) = fp080_flow(j,peak_ti_080(j,1));
        peak_ch080(j,1) = ch_flow080(j,peak_ti_080(j,1));
        fp2ch_frac080(j,1) = peak_fp080(j,1) ./ peak_ch080(j,1);
        fp2tot_frac080(j,1) = peak_fp080(j,1) ./ ...
            (peak_ch080(j,1) + peak_fp080(j,1));
    end
    
    above_base_hold = xs_fp100_cell{i,4}';
    rfp_flow = above_base_hold(xs_rcl_ids(1:324,1)==1,:);
    lfp_flow = above_base_hold(xs_rcl_ids(1:324,1)==3,:);
    fp100_flow = rfp_flow + lfp_flow;
    ch_flow100 = above_base_hold(xs_rcl_ids(1:324,1)==2,:);
    clear above_base_hold rfp_flow lfp_flow
    peak_idtQ_100 = xs_fp100_cell{i,8};
    peak_t_100 = peak_idtQ_100(:,2);
    for j = 1:max_dist/250
        [peak_ti_100(j,1),~] = find(peak_t_100(j,1)==all_time(:,1));
        peak_fp100(j,1) = fp100_flow(j,peak_ti_100(j,1));
        peak_ch100(j,1) = ch_flow100(j,peak_ti_100(j,1));
        fp2ch_frac100(j,1) = peak_fp100(j,1) ./ peak_ch100(j,1);
        fp2tot_frac100(j,1) = peak_fp100(j,1) ./ ...
            (peak_ch100(j,1) + peak_fp100(j,1));
    end
    
%     subplot(3,2,i);
%     scatter(xs_dist,fp2ch_frac040);
%     hold on
%     scatter(xs_dist,fp2ch_frac060);
%     scatter(xs_dist,fp2ch_frac080);
%     scatter(xs_dist,fp2ch_frac100);
%     ylim([0 1.1]);
    
%     subplot(3,2,i);
%     scatter(fp2ch_frac040/max(fp2ch_frac040),fp2ch_frac040);
%     hold on
%     scatter(fp2ch_frac060/max(fp2ch_frac060),fp2ch_frac060);
%     scatter(fp2ch_frac080/max(fp2ch_frac080),fp2ch_frac080);
%     scatter(fp2ch_frac100/max(fp2ch_frac100),fp2ch_frac100);
%     ylim([0 1.2]);
%     xlim([0 1]);

%     subplot(3,2,i);
%     scatter(peak_fp040/max(peak_fp040),fp2ch_frac040);
%     hold on
%     scatter(peak_fp060/max(peak_fp060),fp2ch_frac060);
%     scatter(peak_fp080/max(peak_fp080),fp2ch_frac080);
%     scatter(peak_fp100/max(peak_fp100),fp2ch_frac100);
%     ylim([0 1.2]);
%     xlim([0 1]);

    subplot(3,2,i);
    scatter((peak_ch040+peak_fp040)/peak_index(i),fp2tot_frac040);
    hold on
    scatter((peak_ch060+peak_fp060)/peak_index(i),fp2tot_frac060);
    scatter((peak_ch080+peak_fp080)/peak_index(i),fp2tot_frac080);
    scatter((peak_ch100+peak_fp100)/peak_index(i),fp2tot_frac100);
    ylim([0 1.2]);
    %xlim([0 1]);
    
%     subplot(3,2,i);
%     scatter(peak_ch040+peak_fp040,fp2tot_frac040);
%     hold on
%     scatter(peak_ch060+peak_fp060,fp2tot_frac060);
%     scatter(peak_ch080+peak_fp080,fp2tot_frac080);
%     scatter(peak_ch100+peak_fp100,fp2tot_frac100);
%     ylim([0 1.2]);
    
%     subplot(3,2,i);
%     scatter(xs_dist,fp2tot_frac040);
%     hold on
%     scatter(xs_dist,fp2tot_frac060);
%     scatter(xs_dist,fp2tot_frac080);
%     scatter(xs_dist,fp2tot_frac100);
%     %ylim([0 1.2]);

    
    legend_str = {'FP040','FP060','FP080','FP100'};
    if i == 1
        legend(legend_str,'location','northeast');
    elseif i == 5 
        xlabel('Discharge (m^3/s)');
    elseif i == 6
        xlabel('Discharge (m^3/s)');
    elseif i == 3
        ylabel('Floodplain fraction of flow');
    end
    %xlim([0 28]);
    ylim([0 1]);
    set(gca,'Fontsize',14);

    
end

%% Plot attenuation ratio
figure(2)
set(gcf,'color','white');
xs_dist = (250:250:max_dist)/1000;
for i = 1:6
    
    fp040_attenuation = xs_fp040_cell{i,9};    
    fp060_attenuation = xs_fp060_cell{i,9};
    fp080_attenuation = xs_fp080_cell{i,9};
    fp100_attenuation = xs_fp100_cell{i,9};
    
    max_dist_id = max_dist/250;
    
    subplot(3,2,i);
    plot(xs_dist(1:max_dist_id),fp040_attenuation(1:max_dist_id,1));
    hold on
    plot(xs_dist(1:max_dist_id),fp060_attenuation(1:max_dist_id,1));
    plot(xs_dist(1:max_dist_id),fp080_attenuation(1:max_dist_id,1));
    plot(xs_dist(1:max_dist_id),fp100_attenuation(1:max_dist_id,1));
    if i == 5 
        xlabel('Downstream distance (km)');
        legend('n = 0.040','n = 0.060','n = 0.080','n = 0.100',...
            'Location','southwest');
    elseif i == 6
        xlabel('Downstream distance (km)');
    elseif i == 3
        ylabel('Discharge Attenuation Ratio');
    end
    xlim([0 28]);
    ylim([0.7 1]);
    set(gca,'Fontsize',14);
    
end

%% Plot peak flow fp flow ratios

% Plot 180 cms
rfp_flow = peak180cms040_Q(xs_rcl_ids(1:324,1)==1,:);
lfp_flow = peak180cms040_Q(xs_rcl_ids(1:324,1)==3,:);
fp040_flow180 = rfp_flow + lfp_flow;
chfp040_flow180 = peak180cms040_Q(xs_rcl_ids(1:324,1)==2,:);
clear rfp_flow lfp_flow

rfp_flow = peak180cms060_Q(xs_rcl_ids(1:324,1)==1,:);
lfp_flow = peak180cms060_Q(xs_rcl_ids(1:324,1)==3,:);
fp060_flow180 = rfp_flow + lfp_flow;
chfp060_flow180 = peak180cms060_Q(xs_rcl_ids(1:324,1)==2,:);
clear rfp_flow lfp_flow

rfp_flow = peak180cms080_Q(xs_rcl_ids(1:324,1)==1,:);
lfp_flow = peak180cms080_Q(xs_rcl_ids(1:324,1)==3,:);
fp080_flow180 = rfp_flow + lfp_flow;
chfp080_flow180 = peak180cms080_Q(xs_rcl_ids(1:324,1)==2,:);
clear rfp_flow lfp_flow

rfp_flow = peak180cms100_Q(xs_rcl_ids(1:324,1)==1,:);
lfp_flow = peak180cms100_Q(xs_rcl_ids(1:324,1)==3,:);
fp100_flow180 = rfp_flow + lfp_flow;
chfp100_flow180 = peak180cms100_Q(xs_rcl_ids(1:324,1)==2,:);
clear rfp_flow lfp_flow

figure(3)
scatter(xs_dist,fp040_flow180./(chfp040_flow180+fp040_flow180));
hold on
scatter(xs_dist,fp060_flow180./(chfp060_flow180+fp060_flow180));
scatter(xs_dist,fp080_flow180./(chfp080_flow180+fp080_flow180));
scatter(xs_dist,fp100_flow180./(chfp100_flow180+fp100_flow180));
    
legend_str = {'B040','B060','B080','B100'};
legend(legend_str,'location','northeast');
xlabel('Discharge (m^3/s)');
ylabel('Floodplain fraction of flow');
ylim([0 1]);
set(gca,'Fontsize',14);

% Plot 300 cms
rfp_flow = peak300cms040_Q(xs_rcl_ids(1:324,1)==1,:);
lfp_flow = peak300cms040_Q(xs_rcl_ids(1:324,1)==3,:);
fp040_flow300 = rfp_flow + lfp_flow;
chfp040_flow300 = peak300cms040_Q(xs_rcl_ids(1:324,1)==2,:);
clear rfp_flow lfp_flow

rfp_flow = peak300cms060_Q(xs_rcl_ids(1:324,1)==1,:);
lfp_flow = peak300cms060_Q(xs_rcl_ids(1:324,1)==3,:);
fp060_flow300 = rfp_flow + lfp_flow;
chfp060_flow300 = peak300cms060_Q(xs_rcl_ids(1:324,1)==2,:);
clear rfp_flow lfp_flow

rfp_flow = peak300cms080_Q(xs_rcl_ids(1:324,1)==1,:);
lfp_flow = peak300cms080_Q(xs_rcl_ids(1:324,1)==3,:);
fp080_flow300 = rfp_flow + lfp_flow;
chfp080_flow300 = peak300cms080_Q(xs_rcl_ids(1:324,1)==2,:);
clear rfp_flow lfp_flow

rfp_flow = peak300cms100_Q(xs_rcl_ids(1:324,1)==1,:);
lfp_flow = peak300cms100_Q(xs_rcl_ids(1:324,1)==3,:);
fp100_flow300 = rfp_flow + lfp_flow;
chfp100_flow300 = peak300cms100_Q(xs_rcl_ids(1:324,1)==2,:);
clear rfp_flow lfp_flow

figure(4)
scatter(xs_dist,fp040_flow300./(chfp040_flow300+fp040_flow300));
hold on
scatter(xs_dist,fp060_flow300./(chfp060_flow300+fp060_flow300));
scatter(xs_dist,fp080_flow300./(chfp080_flow300+fp080_flow300));
scatter(xs_dist,fp100_flow300./(chfp100_flow300+fp100_flow300));
    
legend_str = {'B040','B060','B080','B100'};

legend(legend_str,'location','northeast');
xlabel('Discharge (m^3/s)');
ylabel('Floodplain fraction of flow');
ylim([0 1]);
set(gca,'Fontsize',14);