% Extracts water surface elevations from channel and floodplain segments
% under steady flow conditions

clear

%% Input files

load(['/Users/Colin/Documents/MATLAB/DFlow/influences_chap2/'...
    'obs_rcl_ids.mat']);
load(['/Users/Colin/Documents/MATLAB/DFlow/influences_chap2/'...
    'obs_xs_ids.mat']);
load(['/Users/Colin/Documents/MATLAB/DFlow/influences_chap2/'...
    'xs_rcl_ids.mat']);
load(['/Users/Colin/Documents/MATLAB/DFlow/influences_chap2/'...
    'xs_xs_ids.mat']);

bank_peak180cms040_file = ...
    ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3083_180cmsbank040/'...
    'abq10a12bf180cmsbank040rst02_0000_20160104_160000_his.nc'];
bank_peak300cms040_file = ...
    ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3084_300cmsbank040/'...
    'abq10a12bf300cmsbank040rst02_0000_20160103_210000_his.nc'];

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

bank_peak180cms100_file = ...
    ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3087_180cmsbank100/'...
    'abq10a12bf180cmsbank100rst02_0000_20160104_130000_his.nc'];
bank_peak300cms100_file = ...
    ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3088_300cmsbank100/'...
    'abq10a12bf300cmsbank100rst02_0000_20160103_200000_his.nc'];

bank_peak180cms080_file = ...
    ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3085_180cmsbank080/'...
    'abq10a12bf180cmsbank080rst02_0000_20160104_140000_his.nc'];
bank_peak300cms080_file = ...
    ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3086_300cmsbank080/'...
    'abq10a12bf300cmsbank080rst02_0000_20160103_210000_his.nc'];

bank_peak180_files = {bank_peak180cms040_file, bank_peak180cms060_file,...
    bank_peak180cms080_file, bank_peak180cms100_file};
bank_peak300_files = {bank_peak300cms040_file, bank_peak300cms060_file,...
    bank_peak300cms080_file, bank_peak300cms100_file};

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

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

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

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

fp_peak180_files = {fp_peak180cms040_file, fp_peak180cms060_file,...
    fp_peak180cms080_file, fp_peak180cms100_file};
fp_peak300_files = {fp_peak300cms040_file, fp_peak300cms060_file,...
    fp_peak300cms080_file, fp_peak300cms100_file};

peak180cms_chdown_file = ...
    ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12chdown/run10003_180cms/'...
    'abq10a12chdown180cmsrst01_0000_20160102_230000_his.nc'];
peak300cms_chdown_file = ...
    ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12chdown/run10004_300cms/'...
    'abq10a12chdown300cmsrst02_0000_20160103_210000_his.nc'];

peak180cms_chup_file = ...
    ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12chup/run11003_180cms/'...
    'abq10a12upf180cmsrst02_0000_20160104_110000_his.nc'];
peak300cms_chup_file = ...
    ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12chup/run11004_300cms/'...
    'abq10a12upf300cmsrst02_0000_20160103_210000_his.nc'];

ch_peak180_files = {peak180cms_chdown_file, peak180cms_chup_file};
ch_peak300_files = {peak300cms_chdown_file, peak300cms_chup_file};

peak180cms_bankdown_file = ...
    ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12bankdown/run12003_180cms/'...
    'abq10a12bankdownf180rst02_0000_20160104_120000_his.nc'];
peak300cms_bankdown_file = ...
    ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12bankdown/run12004_300cms/'...
    'abq10a12bankdownf300rst02_0000_20160103_200000_his.nc'];

peak180cms_bankup_file = ...
    ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12bankup/run13003_180cms/'...
    'abq10a12bankupf180rst02_0000_20160104_120000_his.nc'];
peak300cms_bankup_file = ...
    ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12bankup/run13004_300cms/'...
    'abq10a12bankupf300rst02_0000_20160103_210000_his.nc'];

banktopo_peak180_files = {peak180cms_bankdown_file, peak180cms_bankup_file};
banktopo_peak300_files = {peak300cms_bankdown_file, peak300cms_bankup_file};

%% Extract data from model runs

num_runs = 4;

for i = 1:num_runs
    
    x_coord = ncread(bank_peak180_files{1,i},'station_x_coordinate');
    y_coord = ncread(bank_peak180_files{1,i},'station_y_coordinate');
    wse_raw = ncread(bank_peak180_files{1,i},'waterlevel');
    wse_end = wse_raw(:,end);
    depth_raw = ncread(bank_peak180_files{1,i},'Waterdepth');
    depth_end = depth_raw(:,end);
    u_raw = ncread(bank_peak180_files{1,i},'x_velocity');
    u_end = u_raw(:,end);
    v_raw = ncread(bank_peak180_files{1,i},'y_velocity');
    v_end = v_raw(:,end);
    Vmag = sqrt(u_end.^2 + v_end.^2);
    bank_raw180_data(i,1) = ...
        {[x_coord y_coord wse_end depth_end u_end v_end Vmag]};
    clear x_coord y_coord wse_raw wse_end depth_raw depth_end...
        u_raw u_end v_raw v_end Vmag
    
    xs_Q_raw = ncread(bank_peak180_files{1,i},'cross_section_discharge');
    xs_Q_bank180(i,1) = {xs_Q_raw(:,end)};
    clear xs_Q_raw
    
    x_coord = ncread(bank_peak300_files{1,i},'station_x_coordinate');
    y_coord = ncread(bank_peak300_files{1,i},'station_y_coordinate');
    wse_raw = ncread(bank_peak300_files{1,i},'waterlevel');
    wse_end = wse_raw(:,end);
    depth_raw = ncread(bank_peak300_files{1,i},'Waterdepth');
    depth_end = depth_raw(:,end);
    u_raw = ncread(bank_peak300_files{1,i},'x_velocity');
    u_end = u_raw(:,end);
    v_raw = ncread(bank_peak300_files{1,i},'y_velocity');
    v_end = v_raw(:,end);
    Vmag = sqrt(u_end.^2 + v_end.^2);
    bank_raw300_data(i,1) = ...
        {[x_coord y_coord wse_end depth_end u_end v_end Vmag]};
    clear x_coord y_coord wse_raw wse_end depth_raw depth_end...
        u_raw u_end v_raw v_end Vmag
    
    xs_Q_raw = ncread(bank_peak300_files{1,i},'cross_section_discharge');
    xs_Q_bank300(i,1) = {xs_Q_raw(:,end)};
    clear xs_Q_raw
    
    x_coord = ncread(fp_peak180_files{1,i},'station_x_coordinate');
    y_coord = ncread(fp_peak180_files{1,i},'station_y_coordinate');
    wse_raw = ncread(fp_peak180_files{1,i},'waterlevel');
    wse_end = wse_raw(:,end);
    depth_raw = ncread(fp_peak180_files{1,i},'Waterdepth');
    depth_end = depth_raw(:,end);
    u_raw = ncread(fp_peak180_files{1,i},'x_velocity');
    u_end = u_raw(:,end);
    v_raw = ncread(fp_peak180_files{1,i},'y_velocity');
    v_end = v_raw(:,end);
    Vmag = sqrt(u_end.^2 + v_end.^2);
    fp_raw180_data(i,1) = ...
        {[x_coord y_coord wse_end depth_end u_end v_end Vmag]};
    clear x_coord y_coord wse_raw wse_end depth_raw depth_end...
        u_raw u_end v_raw v_end Vmag
    
    xs_Q_raw = ncread(fp_peak180_files{1,i},'cross_section_discharge');
    xs_Q_fp180(i,1) = {xs_Q_raw(:,end)};
    clear xs_Q_raw
    
    x_coord = ncread(fp_peak300_files{1,i},'station_x_coordinate');
    y_coord = ncread(fp_peak300_files{1,i},'station_y_coordinate');
    wse_raw = ncread(fp_peak300_files{1,i},'waterlevel');
    wse_end = wse_raw(:,end);
    depth_raw = ncread(fp_peak300_files{1,i},'Waterdepth');
    depth_end = depth_raw(:,end);
    u_raw = ncread(fp_peak300_files{1,i},'x_velocity');
    u_end = u_raw(:,end);
    v_raw = ncread(fp_peak300_files{1,i},'y_velocity');
    v_end = v_raw(:,end);
    Vmag = sqrt(u_end.^2 + v_end.^2);
    fp_raw300_data(i,1) = ...
        {[x_coord y_coord wse_end depth_end u_end v_end Vmag]};
    clear x_coord y_coord wse_raw wse_end depth_raw depth_end...
        u_raw u_end v_raw v_end Vmag
    
    xs_Q_raw = ncread(fp_peak300_files{1,i},'cross_section_discharge');
    xs_Q_fp300(i,1) = {xs_Q_raw(:,end)};
    clear xs_Q_raw
    
end

for i = 1:2
    
    x_coord = ncread(ch_peak180_files{1,i},'station_x_coordinate');
    y_coord = ncread(ch_peak180_files{1,i},'station_y_coordinate');
    wse_raw = ncread(ch_peak180_files{1,i},'waterlevel');
    wse_end = wse_raw(:,end);
    depth_raw = ncread(ch_peak180_files{1,i},'Waterdepth');
    depth_end = depth_raw(:,end);
    u_raw = ncread(ch_peak180_files{1,i},'x_velocity');
    u_end = u_raw(:,end);
    v_raw = ncread(ch_peak180_files{1,i},'y_velocity');
    v_end = v_raw(:,end);
    Vmag = sqrt(u_end.^2 + v_end.^2);
    ch_raw180_data(i,1) = ...
        {[x_coord y_coord wse_end depth_end u_end v_end Vmag]};
    clear x_coord y_coord wse_raw wse_end depth_raw depth_end...
        u_raw u_end v_raw v_end Vmag
    
    xs_Q_raw = ncread(ch_peak180_files{1,i},'cross_section_discharge');
    xs_Q_ch180(i,1) = {xs_Q_raw(:,end)};
    clear xs_Q_raw
    
    x_coord = ncread(ch_peak300_files{1,i},'station_x_coordinate');
    y_coord = ncread(ch_peak300_files{1,i},'station_y_coordinate');
    wse_raw = ncread(ch_peak300_files{1,i},'waterlevel');
    wse_end = wse_raw(:,end);
    depth_raw = ncread(ch_peak300_files{1,i},'Waterdepth');
    depth_end = depth_raw(:,end);
    u_raw = ncread(ch_peak300_files{1,i},'x_velocity');
    u_end = u_raw(:,end);
    v_raw = ncread(ch_peak300_files{1,i},'y_velocity');
    v_end = v_raw(:,end);
    Vmag = sqrt(u_end.^2 + v_end.^2);
    ch_raw300_data(i,1) = ...
        {[x_coord y_coord wse_end depth_end u_end v_end Vmag]};
    clear x_coord y_coord wse_raw wse_end depth_raw depth_end...
        u_raw u_end v_raw v_end Vmag
    
    xs_Q_raw = ncread(ch_peak300_files{1,i},'cross_section_discharge');
    xs_Q_ch300(i,1) = {xs_Q_raw(:,end)};
    clear xs_Q_raw
  
    x_coord = ncread(banktopo_peak180_files{1,i},'station_x_coordinate');
    y_coord = ncread(banktopo_peak180_files{1,i},'station_y_coordinate');
    wse_raw = ncread(banktopo_peak180_files{1,i},'waterlevel');
    wse_end = wse_raw(:,end);
    depth_raw = ncread(banktopo_peak180_files{1,i},'Waterdepth');
    depth_end = depth_raw(:,end);
    u_raw = ncread(banktopo_peak180_files{1,i},'x_velocity');
    u_end = u_raw(:,end);
    v_raw = ncread(banktopo_peak180_files{1,i},'y_velocity');
    v_end = v_raw(:,end);
    Vmag = sqrt(u_end.^2 + v_end.^2);
    banktopo_raw180_data(i,1) = ...
        {[x_coord y_coord wse_end depth_end u_end v_end Vmag]};
    clear x_coord y_coord wse_raw wse_end depth_raw depth_end...
        u_raw u_end v_raw v_end Vmag
    
    xs_Q_raw = ncread(banktopo_peak180_files{1,i},'cross_section_discharge');
    xs_Q_banktopo180(i,1) = {xs_Q_raw(:,end)};
    clear xs_Q_raw
    
    x_coord = ncread(banktopo_peak300_files{1,i},'station_x_coordinate');
    y_coord = ncread(banktopo_peak300_files{1,i},'station_y_coordinate');
    wse_raw = ncread(banktopo_peak300_files{1,i},'waterlevel');
    wse_end = wse_raw(:,end);
    depth_raw = ncread(banktopo_peak300_files{1,i},'Waterdepth');
    depth_end = depth_raw(:,end);
    u_raw = ncread(banktopo_peak300_files{1,i},'x_velocity');
    u_end = u_raw(:,end);
    v_raw = ncread(banktopo_peak300_files{1,i},'y_velocity');
    v_end = v_raw(:,end);
    Vmag = sqrt(u_end.^2 + v_end.^2);
    banktopo_raw300_data(i,1) = ...
        {[x_coord y_coord wse_end depth_end u_end v_end Vmag]};
    clear x_coord y_coord wse_raw wse_end depth_raw depth_end...
        u_raw u_end v_raw v_end Vmag
    
    xs_Q_raw = ncread(...
        banktopo_peak300_files{1,i},'cross_section_discharge');
    xs_Q_banktopo300(i,1) = {xs_Q_raw(:,end)};
    clear xs_Q_raw
    
end


%% Average Water Surface Elevation in Channel

max_xs = 108;
bank_wse180 = zeros(max_xs,4);
bank_wse300 = zeros(max_xs,4);
bank_velocity180 = zeros(max_xs,4);
bank_velocity300 = zeros(max_xs,4);

fp_wse180 = zeros(max_xs,4);
fp_wse300 = zeros(max_xs,4);
fp_velocity180 = zeros(max_xs,4);
fp_velocity300 = zeros(max_xs,4);

prop_fp_ch_180fp = zeros(max_xs,4);
prop_fp_ch_300fp = zeros(max_xs,4);
prop_fp_ch_180bank = zeros(max_xs,4);
prop_fp_ch_300bank = zeros(max_xs,4);

for i = 1:4
    
    bank_obs_data180 = bank_raw180_data{i,1};
    bank_obs_data300 = bank_raw300_data{i,1};
    
    fp_obs_data180 = fp_raw180_data{i,1};
    fp_obs_data300 = fp_raw300_data{i,1};
    
    bank_xs_data180 = xs_Q_bank180{i,1};
    bank_xs_data300 = xs_Q_bank300{i,1};
    
    fp_xs_data180 = xs_Q_fp180{i,1};
    fp_xs_data300 = xs_Q_fp300{i,1};
    
    for j = 1:max_xs
        
        bank_wse180(j,i) = mean(bank_obs_data180(obs_rcl_ids==2 & ...
            obs_xs_ids==j & bank_obs_data180(:,4)>0,3));
        bank_wse300(j,i) = mean(bank_obs_data300(obs_rcl_ids==2 &...
            obs_xs_ids==j & bank_obs_data300(:,4)>0,3));
        bank_velocity180(j,i) = mean(bank_obs_data180(obs_rcl_ids==2&...
            obs_xs_ids==j&bank_obs_data180(:,4)>0,7));
        bank_velocity300(j,i) = mean(bank_obs_data300(obs_rcl_ids==2&...
            obs_xs_ids==j&bank_obs_data300(:,4)>0,7));
        
        fp_wse180(j,i) = mean(fp_obs_data180(obs_rcl_ids==2&...
            obs_xs_ids==j & fp_obs_data180(:,4)>0,3));
        fp_wse300(j,i) = mean(fp_obs_data300(obs_rcl_ids==2&...
            obs_xs_ids==j & fp_obs_data300(:,4)>0,3));
        fp_velocity180(j,i) = mean(fp_obs_data180(obs_rcl_ids==2&...
            obs_xs_ids==j&fp_obs_data180(:,4)>0,7));
        fp_velocity300(j,i) = mean(fp_obs_data300(obs_rcl_ids==2&...
            obs_xs_ids==j&fp_obs_data300(:,4)>0,7));
        
        prop_fp_ch_180fp(j,i) = sum(fp_xs_data180(xs_rcl_ids~=2 & ...
            xs_xs_ids==j)) / sum(fp_xs_data180(xs_xs_ids==j));
        prop_fp_ch_300fp(j,i) = sum(fp_xs_data300(xs_rcl_ids~=2 & ...
            xs_xs_ids==j)) / sum(fp_xs_data300(xs_xs_ids==j));
        
        prop_fp_ch_180bank(j,i) = sum(bank_xs_data180(xs_rcl_ids~=2 & ...
            xs_xs_ids==j)) / sum(bank_xs_data180(xs_xs_ids==j));
        prop_fp_ch_300bank(j,i) = sum(bank_xs_data300(xs_rcl_ids~=2 & ...
            xs_xs_ids==j)) / sum(bank_xs_data300(xs_xs_ids==j));
        
    end
    
end

ch_wse180 = zeros(max_xs,2);
ch_wse300 = zeros(max_xs,2);
ch_velocity180 = zeros(max_xs,2);
ch_velocity300 = zeros(max_xs,2);
prop_fp_ch_180ch = zeros(max_xs,2);
prop_fp_ch_300ch = zeros(max_xs,2);

banktopo_wse180 = zeros(max_xs,2);
banktopo_wse300 = zeros(max_xs,2);
banktopo_velocity180 = zeros(max_xs,2);
banktopo_velocity300 = zeros(max_xs,2);
prop_fp_ch_180banktopo = zeros(max_xs,2);
prop_fp_ch_300banktopo = zeros(max_xs,2);

for i = 1:2

    ch_obs_data180 = ch_raw180_data{i,1};
    ch_obs_data300 = ch_raw300_data{i,1};
    
    ch_xs_data180 = xs_Q_ch180{i,1};
    ch_xs_data300 = xs_Q_ch300{i,1};
    
    banktopo_obs_data180 = banktopo_raw180_data{i,1};
    banktopo_obs_data300 = banktopo_raw300_data{i,1};
    
    banktopo_xs_data180 = xs_Q_banktopo180{i,1};
    banktopo_xs_data300 = xs_Q_banktopo300{i,1};
    
    for j = 1:max_xs
        
        ch_wse180(j,i) = mean(ch_obs_data180(obs_rcl_ids==2&...
            obs_xs_ids==j & ch_obs_data180(:,4)>0,3));
        ch_wse300(j,i) = mean(ch_obs_data300(obs_rcl_ids==2&...
            obs_xs_ids==j & ch_obs_data300(:,4)>0,3));
        ch_velocity180(j,i) = mean(ch_obs_data180(obs_rcl_ids==2&...
            obs_xs_ids==j&ch_obs_data180(:,4)>0,7));
        ch_velocity300(j,i) = mean(ch_obs_data300(obs_rcl_ids==2&...
            obs_xs_ids==j&ch_obs_data300(:,4)>0,7));
        
        prop_fp_ch_180ch(j,i) = sum(ch_xs_data180(xs_rcl_ids~=2 & ...
            xs_xs_ids==j)) / sum(ch_xs_data180(xs_xs_ids==j));
        prop_fp_ch_300ch(j,i) = sum(ch_xs_data300(xs_rcl_ids~=2 & ...
            xs_xs_ids==j)) / sum(ch_xs_data300(xs_xs_ids==j));
        
        banktopo_wse180(j,i) = mean(...
            banktopo_obs_data180(obs_rcl_ids==2&...
            obs_xs_ids==j & banktopo_obs_data180(:,4)>0,3));
        banktopo_wse300(j,i) = mean(...
            banktopo_obs_data300(obs_rcl_ids==2&...
            obs_xs_ids==j & banktopo_obs_data300(:,4)>0,3));
        banktopo_velocity180(j,i) = mean(...
            banktopo_obs_data180(obs_rcl_ids==2&...
            obs_xs_ids==j&banktopo_obs_data180(:,4)>0,7));
        banktopo_velocity300(j,i) = mean(...
            banktopo_obs_data300(obs_rcl_ids==2&...
            obs_xs_ids==j&banktopo_obs_data300(:,4)>0,7));
        
        prop_fp_ch_180banktopo(j,i) = sum(...
            banktopo_xs_data180(xs_rcl_ids~=2 & ...
            xs_xs_ids==j)) / sum(banktopo_xs_data180(xs_xs_ids==j));
        prop_fp_ch_300banktopo(j,i) = sum(...
            banktopo_xs_data300(xs_rcl_ids~=2 & ...
            xs_xs_ids==j)) / sum(banktopo_xs_data300(xs_xs_ids==j));
        
    end
    
end

%% Average Water Surface Elevation on FP
% 
% max_xs = 108;
% fp_bank_wse180 = zeros(max_xs,4);
% fp_bank_wse300 = zeros(max_xs,4);
% fp_bank_velocity180 = zeros(max_xs,4);
% fp_bank_velocity300 = zeros(max_xs,4);
% 
% fp_fp_wse180 = zeros(max_xs,4);
% fp_fp_wse300 = zeros(max_xs,4);
% fp_fp_velocity180 = zeros(max_xs,4);
% fp_fp_velocity300 = zeros(max_xs,4);
% 
% for i = 1:4
%     
%     bank_obs_data180 = bank_raw180_data{i,1};
%     bank_obs_data300 = bank_raw300_data{i,1};
%     
%     fp_obs_data180 = fp_raw180_data{i,1};
%     fp_obs_data300 = fp_raw300_data{i,1};
%     
%     for j = 1:max_xs
%         
%         fp_bank_wse180(j,i) = mean(bank_obs_data180((...
%             obs_rcl_ids==1|obs_rcl_ids==3)&...
%             obs_xs_ids==j & bank_obs_data180(:,4)>0,3));
%         fp_bank_wse300(j,i) = mean(bank_obs_data300((...
%             obs_rcl_ids==1|obs_rcl_ids==3)&...
%             obs_xs_ids==j & bank_obs_data300(:,4)>0,3));
%         fp_bank_velocity180(j,i) = mean(bank_obs_data180((...
%             obs_rcl_ids==1|obs_rcl_ids==3)&...
%             obs_xs_ids==j&bank_obs_data180(:,4)>0,7));
%         fp_bank_velocity300(j,i) = mean(bank_obs_data300((...
%             obs_rcl_ids==1|obs_rcl_ids==3)&...
%             obs_xs_ids==j&bank_obs_data300(:,4)>0,7));
%         
%         fp_fp_wse180(j,i) = mean(fp_obs_data180((...
%             obs_rcl_ids==1|obs_rcl_ids==3)&...
%             obs_xs_ids==j & fp_obs_data180(:,4)>0,3));
%         fp_fp_wse300(j,i) = mean(fp_obs_data300((...
%             obs_rcl_ids==1|obs_rcl_ids==3)&...
%             obs_xs_ids==j & fp_obs_data300(:,4)>0,3));
%         fp_fp_velocity180(j,i) = mean(fp_obs_data180((...
%             obs_rcl_ids==1|obs_rcl_ids==3)&...
%             obs_xs_ids==j&fp_obs_data180(:,4)>0,7));
%         fp_fp_velocity300(j,i) = mean(fp_obs_data300((...
%             obs_rcl_ids==1|obs_rcl_ids==3)&...
%             obs_xs_ids==j&fp_obs_data300(:,4)>0,7));
%         
%     end
%     
% end

%% Plot mean channel elevation difference from Baseline

a = 0;
for i = [1 3 4]
    a = a + 1;
    diff_fp_wse180(:,a) = fp_wse180(:,i) - fp_wse180(:,2);
    diff_fp_wse300(:,a) = fp_wse300(:,i) - fp_wse300(:,2);
    diff_bank_wse180(:,a) = bank_wse180(:,i) - fp_wse180(:,2);
    diff_bank_wse300(:,a) = bank_wse300(:,i) - fp_wse300(:,2);
    
    diff_fp_velocity180(:,a) = fp_velocity180(:,i) - fp_velocity180(:,2);
    diff_fp_velocity300(:,a) = fp_velocity300(:,i) - fp_velocity300(:,2);
    diff_bank_velocity180(:,a) = bank_velocity180(:,i) - ...
        fp_velocity180(:,2);
    diff_bank_velocity300(:,a) = bank_velocity300(:,i) - ...
        fp_velocity300(:,2);
    
    diff_prop_fp180(:,a) = prop_fp_ch_180fp(:,i) - prop_fp_ch_180fp(:,2);
    diff_prop_fp300(:,a) = prop_fp_ch_300fp(:,i) - prop_fp_ch_300fp(:,2);
    diff_prop_bank180(:,a) = prop_fp_ch_180bank(:,i) - ...
        prop_fp_ch_180fp(:,2);
    diff_prop_bank300(:,a) = prop_fp_ch_300bank(:,i) - ...
        prop_fp_ch_300fp(:,2);
    
%     diff_fp_fp_wse180(:,a) = fp_fp_wse180(:,i) - fp_fp_wse180(:,2);
%     diff_fp_fp_wse300(:,a) = fp_fp_wse300(:,i) - fp_fp_wse300(:,2);
    if i == 1 || i == 3
        diff_ch_wse180(:,a) = ch_wse180(:,a) - fp_wse180(:,2);
        diff_ch_wse300(:,a) = ch_wse300(:,a) - fp_wse300(:,2);
        diff_banktopo_wse180(:,a) = banktopo_wse180(:,a) - fp_wse180(:,2);
        diff_banktopo_wse300(:,a) = banktopo_wse300(:,a) - fp_wse300(:,2);
        
        diff_ch_velocity180(:,a) = ch_velocity180(:,a) -...
            fp_velocity180(:,2);
        diff_ch_velocity300(:,a) = ch_velocity300(:,a) -... 
            fp_velocity300(:,2);
        diff_banktopo_velocity180(:,a) = banktopo_velocity180(:,a) -...
            fp_velocity180(:,2);
        diff_banktopo_velocity300(:,a) = banktopo_velocity300(:,a) -... 
            fp_velocity300(:,2);
        
        diff_prop_ch180(:,a) = prop_fp_ch_180ch(:,a) - ...
            prop_fp_ch_180fp(:,2);
        diff_prop_ch300(:,a) = prop_fp_ch_300ch(:,a) - ...
            prop_fp_ch_300fp(:,2);
        diff_prop_banktopo180(:,a) = prop_fp_ch_180banktopo(:,a) - ...
            prop_fp_ch_180fp(:,2);
        diff_prop_banktopo300(:,a) = prop_fp_ch_300banktopo(:,a) - ...
            prop_fp_ch_300fp(:,2);
    end
end

mean_fp100_wse = mean(diff_fp_wse300(:,3));
mean_chup_wse = mean(diff_ch_wse300(:,2));
mean_fp040_wse = mean(diff_fp_wse300(:,1));
mean_fp100_V180 = mean(diff_fp_velocity180(:,3));
mean_fp100_V300 = mean(diff_fp_velocity300(:,3));
mean_fp040_V180 = mean(diff_fp_velocity180(:,1));
mean_fp040_V300 = mean(diff_fp_velocity300(:,1));
mean_V180_n = mean_fp040_V180 - mean_fp100_V180;
mean_V300_n = mean_fp040_V300 - mean_fp100_V300;
mean_prop_bank040_300 = mean(diff_prop_bank300(:,1));
mean_prop_bank100_300 = mean(diff_prop_bank300(:,3));
mean_prop_bank = mean_prop_bank040_300 - mean_prop_bank100_300;
mean_prop_fp040_300 = mean(diff_prop_fp300(:,1));
mean_prop_fp100_300 = mean(diff_prop_fp300(:,3));
mean_prop_fp = mean_prop_fp040_300 - mean_prop_fp100_300;

%% Steady Subplots

h1 = figure(1);
h1.Units = 'inches';
h1.PaperSize = [8.5 11];
set(gcf,'color','white');
%label_str3 = {'CH Down','CH Up','FP040','FP080','FP100'};
label_str3 = {'','','','','','','','','',''};
spheight = 0.13;
gap_size = 0.05;

sp1=subplot(6,1,1);
p1 = get(sp1, 'pos');
boxplot([diff_ch_wse180 diff_banktopo_wse180 diff_fp_wse180 ...
    diff_bank_wse180],...
    'Labels',label_str3,'Symbol','k.');
p1(2) = 0.1+spheight*5+gap_size*2;
p1(4) = spheight;
set(sp1, 'pos', p1);
set(gca,'Fontsize',14);
% ylabel('Relative Depth (m)');

sp2=subplot(6,1,2);
p2 = get(sp2, 'pos');
boxplot([diff_ch_wse300 diff_banktopo_wse300...
    diff_fp_wse300  diff_bank_wse300],...
    'Labels',label_str3,'Symbol','k.');
p2(2) = 0.1+spheight*4+gap_size*2;
p2(4) = spheight;
set(sp2, 'pos', p2);
set(gca,'Fontsize',14);
ylabel({'Relative WSE (m)';...
    '  Major           Moderate'});
ylab = get(gca,'YLabel');
set(ylab,'Position',get(ylab,'Position') - [0 -0.3 0])

sp3=subplot(6,1,3);
p3 = get(sp3, 'pos');
boxplot([diff_ch_velocity180 diff_banktopo_velocity180...
    diff_fp_velocity180 diff_bank_velocity180],...
    'Labels',label_str3,'Symbol','k.')
p3(2) = 0.1+spheight*3+gap_size;
p3(4) = spheight;
set(sp3, 'pos', p3);
set(gca,'Fontsize',14);
ylim([-0.3 0.3]);
yticks(-0.2:0.2:0.2);
% ylabel('Relative Velocity (m/s)');

sp4=subplot(6,1,4);
p4 = get(sp4, 'pos');
boxplot([diff_ch_velocity300 diff_banktopo_velocity300...
    diff_fp_velocity300 diff_bank_velocity300],...
    'Labels',label_str3,'Symbol','k.')
p4(2) = 0.1+spheight*2+gap_size;
p4(4) = spheight;
set(sp4, 'pos', p4);
set(gca,'Fontsize',14);
ylim([-0.3 0.3]);
yticks(-0.2:0.2:0.2);
ylabel({'Relative Velocity (ms^{ -1})';'  Major           Moderate'});
ylab = get(gca,'YLabel');
set(ylab,'Position',get(ylab,'Position') - [0 -0.3 0])

sp5=subplot(6,1,5);
p5 = get(sp5, 'pos');
boxplot([diff_prop_ch180 diff_prop_banktopo180...
    diff_prop_fp180 diff_prop_bank180],...
    'Labels',label_str3,'Symbol','k.')
p5(2) = 0.1+spheight;
p5(4) = spheight;
set(sp5, 'pos', p5);
set(gca,'Fontsize',14);
% ylabel('Relative Q_{FP}/Q_{Total}');

sp6=subplot(6,1,6);
p6 = get(sp6, 'pos');
boxplot([diff_prop_ch300 diff_prop_banktopo300...
    diff_prop_fp300 diff_prop_bank300],...
    'Labels',label_str3,'Symbol','k.')
p6(2) = 0.1;
p6(4) = spheight;
set(sp6, 'pos', p6);
set(gca,'Fontsize',14);
ylabel({'Relative Q_{FP}/Q_{Total}';'  Major          Moderate'});
ylab = get(gca,'YLabel');
set(ylab,'Position',get(ylab,'Position') - [0 -0.14 0])

print('-fillpage',['/Users/Colin/Documents/MATLAB/DFlow/'...
    'influences_chap2/figures/steady_subplots'],'-dpdf');

%%

% h2 = figure(2);
% h2.Units = 'inches';
% h2.PaperSize = [5 5.5];
% set(gcf,'color','white');
% % label_str3 = {'CH Down','CH Up','FP040','FP080','FP100',...
% %     'Bank040','Bank080','Bank100'};
% label_str3 = {'','','','','','','',''};
% subplot(2,1,1)
% boxplot(...
%     [diff_ch_wse180 diff_fp_wse180 diff_bank_wse180],...
%     'Labels',label_str3,'Symbol','k.')
% set(gca,'Fontsize',14);
% ylabel('Relative Depth (m)');
% 
% subplot(2,1,2)
% boxplot(...
%     [diff_ch_wse300 diff_fp_wse300 diff_bank_wse300],...
%     'Labels',label_str3,'Symbol','k.')
% set(gca,'Fontsize',14);
% ylabel('Relative Depth (m)');
% 
% print('-fillpage',['/Users/Colin/Documents/MATLAB/DFlow/'...
%     'influences_chap2/figures/ch_depth'],'-dpdf');
% 
% % Differences in velocity 
% h3 = figure(3);
% h3.Units = 'inches';
% h3.PaperSize = [5 5.5];
% set(gcf,'color','white');
% % label_str3 = {'CH Down','CH Up','FP040','FP080','FP100',...
% %     'Bank040','Bank080','Bank100'};
% label_str3 = {'','','','','','','',''};
% subplot(2,1,1)
% boxplot([diff_ch_velocity180 diff_fp_velocity180 diff_bank_velocity180],...
%     'Labels',label_str3,'Symbol','k.')
% set(gca,'Fontsize',14);
% ylabel('Relative Velocity (m/s)');
% 
% subplot(2,1,2)
% boxplot([diff_ch_velocity300 diff_fp_velocity300 diff_bank_velocity300],...
%     'Labels',label_str3,'Symbol','k.')
% set(gca,'Fontsize',14);
% ylabel('Relative Velocity (m/s)');
% 
% print('-fillpage',['/Users/Colin/Documents/MATLAB/DFlow/'...
%     'influences_chap2/figures/ch_velocity'],'-dpdf');
% 
% % Proportion of floodplain flow to total flow
% h4 = figure(4);
% h4.Units = 'inches';
% h4.PaperSize = [5 5.5];
% set(gcf,'color','white');
% label_str3 = {'','','','','','','',''};
% subplot(2,1,1)
% boxplot([diff_prop_ch180 diff_prop_fp180 diff_prop_bank180],...
%     'Labels',label_str3,'Symbol','k.')
% set(gca,'Fontsize',14);
% ylabel('Relative Q_{FP}/Q_{Total}');
% 
% subplot(2,1,2)
% boxplot([diff_prop_ch300 diff_prop_fp300 diff_prop_bank300],...
%     'Labels',label_str3,'Symbol','k.')
% set(gca,'Fontsize',14);
% ylabel('Relative Q_{FP}/Q_{Total}');
% 
% print('-fillpage',['/Users/Colin/Documents/MATLAB/DFlow/'...
%     'influences_chap2/figures/proportion_Qfp'],'-dpdf');

% subplot(2,2,2)
% label_str = {'FP040','FP080','FP100'};
% label_str2 = {'CH Down','CH Up'};
% boxplot(diff_fp_wse180,'Labels',label_str)
% ylim([-0.2 0.3]);
% 
% subplot(2,2,1)
% boxplot(diff_ch_wse180,'Labels',label_str2)
% %ylim([-0.2 0.3]);
% 
% subplot(2,2,4)
% boxplot(diff_fp_wse300,'Labels',label_str)
% ylim([-0.2 0.3]);
% 
% subplot(2,2,3)
% boxplot(diff_ch_wse300,'Labels',label_str2);
% %ylim([-0.2 0.3]);

%% Plot extras

% down_dist = (250:250:27000)/1000;
% 
% figure(1)
% set(gcf,'color','white');
% plot(down_dist,bank_wse180(:,1));
% hold on
% plot(down_dist,bank_wse180(:,2));
% plot(down_dist,bank_wse180(:,3));
% plot(down_dist,bank_wse180(:,4));
% ax = gca;
% ax.ColorOrderIndex = 1;
% plot(down_dist,fp_wse180(:,1),'--');
% plot(down_dist,fp_wse180(:,2),'--');
% plot(down_dist,fp_wse180(:,3),'--');
% plot(down_dist,fp_wse180(:,4),'--');
% legend('0.040', '0.060','0.080','0.100')
% xlim([10 15]);
% 
% figure(2)
% set(gcf,'color','white');
% plot(down_dist,bank_velocity180(:,1));
% hold on
% plot(down_dist,bank_velocity180(:,2));
% plot(down_dist,bank_velocity180(:,3));
% plot(down_dist,bank_velocity180(:,4));
% ax = gca;
% ax.ColorOrderIndex = 1;
% plot(down_dist,bank_velocity300(:,1),'--');
% plot(down_dist,bank_velocity300(:,2),'--');
% plot(down_dist,bank_velocity300(:,3),'--');
% plot(down_dist,bank_velocity300(:,4),'--');
% xlim([10 15]);
% 
% figure(3)
% set(gcf,'color','white');
% plot(down_dist,fp_wse180(:,1));
% hold on
% plot(down_dist,fp_wse180(:,2));
% plot(down_dist,fp_wse180(:,3));
% plot(down_dist,fp_wse180(:,4));
% ax = gca;
% ax.ColorOrderIndex = 1;
% plot(down_dist,fp_wse300(:,1),'--');
% plot(down_dist,fp_wse300(:,2),'--');
% plot(down_dist,fp_wse300(:,3),'--');
% plot(down_dist,fp_wse300(:,4),'--');
% xlim([5 15]);
% 
% figure(4)
% set(gcf,'color','white');
% plot(down_dist,fp_velocity180(:,1));
% hold on
% plot(down_dist,fp_velocity180(:,2));
% plot(down_dist,fp_velocity180(:,3));
% plot(down_dist,fp_velocity180(:,4));
% ax = gca;
% ax.ColorOrderIndex = 1;
% plot(down_dist,fp_velocity300(:,1),'--');
% plot(down_dist,fp_velocity300(:,2),'--');
% plot(down_dist,fp_velocity300(:,3),'--');
% plot(down_dist,fp_velocity300(:,4),'--');
% xlim([5 15]);

% figure(5)
% bank_fp_Vmag180 = [bank_velocity180(:,1) fp_velocity180(:,1)...
%     bank_velocity180(:,2) fp_velocity180(:,2)...
%     bank_velocity180(:,3) fp_velocity180(:,3)...
%     bank_velocity180(:,4) fp_velocity180(:,4)];
% bank_fp_Vmag300 = [bank_velocity300(:,1) fp_velocity300(:,1)...
%     bank_velocity300(:,2) fp_velocity300(:,2)...
%     bank_velocity300(:,3) fp_velocity300(:,3)...
%     bank_velocity300(:,4) fp_velocity300(:,4)];
% subplot(2,1,1)
% label_str = {'B040','FP040','B060','FP060','B080','FP080','B100','FP100'};
% boxplot(bank_fp_Vmag180,'PlotStyle','compact',...
%     'Labels',label_str)
% 
% subplot(2,1,2)
% boxplot(bank_fp_Vmag300,'PlotStyle','compact',...
%     'Labels',label_str)