% Evaluate peak water surface elevations for distributed roughness model in
% comparison with FP/Bank elevation

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];
max_xs_num = max_dist/250;

% Load observation ids
id_floc = '/Users/Colin/Documents/MATLAB/DFlow/influences_chap2/';
load(strcat(id_floc,'obs_xs_ids.mat'));
load(strcat(id_floc,'obs_rcl_ids.mat'));

% Load xs minimum elevations
load(strcat(id_floc,'data/xs_elevation.mat'));
xs_elevation = xs_elevation(2:end,:);

%% 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_wse = ncread(base30cms040_file,'waterlevel');
% base150cms040_all_wse = ncread(base150cms040_file,'waterlevel');
% 
% % Extract steady state (end time) flow at all cross-sections
% base30cms040_wse = base30cms040_all_wse(:,end);
% base150cms040_wse = base150cms040_all_wse(:,end);
% 
% clear base30cms040_all_wse base150cms040_all_wse

%% Collect flow under steady peak flow for distributed n FP
% 
% peak180cms_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
%     'abq2010add2012bath/run3045_180cmsdist/'...
%     'abq10a12bf180cmsrst02_0000_20160104_110000_his.nc'];
% peak300cms_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
%     'abq2010add2012bath/run3046_300cmsdist/'...
%     'abq10a12bf300cmsrst02_0000_20160103_210000_his.nc'];
% 
% % Extract WSE at all cross-section times
% peak180cms_all_wse = ncread(peak180cms_file,'waterlevel');
% peak300cms_all_wse = ncread(peak300cms_file,'waterlevel');
% 
% % Extract steady state (end time) flow at all cross-sections
% peak180cms_wse_steady = peak180cms_all_wse(:,end);
% peak300cms_wse_steady = peak300cms_all_wse(:,end);
% 
% % Extract WSE at all cross-section times
% peak180cms_all_depth = ncread(peak180cms_file,'Waterdepth');
% peak300cms_all_depth = ncread(peak300cms_file,'Waterdepth');
% 
% % Extract steady state (end time) flow at all cross-sections
% peak180cms_depth_steady = peak180cms_all_depth(:,end);
% peak300cms_depth_steady = peak300cms_all_depth(:,end);
% 
% clear peak180cms_all_wse peak300cms_all_wse ...
%     peak180cms_all_depth peak300cms_all_depth
% 
% % Find maximum WSE for peak flow at each XS
% peak180cms_wse = zeros(max_xs_num,1);
% peak300cms_wse = zeros(max_xs_num,1);
% for i = 1:max_xs_num
%     peak180cms_wse(i,1) = mean(...
%         peak180cms_wse_steady(obs_rcl_ids(:,1)==2 &...
%         obs_xs_ids(:,1)==i & peak180cms_depth_steady(:,1)>0,1));
%     peak300cms_wse(i,1) = mean(...
%         peak300cms_wse_steady(obs_rcl_ids(:,1)==2 &...
%         obs_xs_ids(:,1)==i & peak300cms_depth_steady(:,1)>0,1));
% end
% 
% %% 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 = ['run3059_s180dist/';'run3060_s300dist/';...
%     'run3061_m180dist/';'run3062_m300dist/';...
%     'run3063_l180dist/';'run3064_l300dist/']; % file path
% 
% folder_num = numel(fpath(:,1));
% obs_fp_cell = cell(folder_num,7);
% 
% for fi = 1:folder_num
%     
%     obs_fp_cell(fi,1) = {fpath(fi,:)};
%     upstream_peak = peak_index(fi);
% 
%     if mod(fi,2) == 0
%         obs_peak_wse = peak300cms_wse;
%         %obs_base_wse = base150cms040_wse;
%     elseif mod(fi,2) == 1
%         obs_peak_wse = peak180cms_wse;
%         %obs_base_wse = base30cms040_wse;
%     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
%     
%     [obs_fp_cell(fi,2:7)] = ...
%         hpeak_obs_from_his_chap2(fnum, fnames, fpath_fi, max_dist,min_t,...
%         max_t, obs_peak_wse, xs_elevation, obs_xs_ids, obs_rcl_ids);
%     
% end

% Load cell array that is computed in commented out section above
load(['/Users/Colin/Documents/MATLAB/DFlow/influences_chap2/SAR_cells/'...
    'obs_fp_cell.mat']);

%% Compare to FP and bank elevation

% Load bank and floodplain elevation
load(['/Users/Colin/Documents/MATLAB/DFlow/influences_chap2/data/'...
    'bank_elevation_1000m.mat']);
load(['/Users/Colin/Documents/MATLAB/DFlow/influences_chap2/data/'...
    'fp_elevation_1000m.mat']);
load(['/Users/Colin/Documents/MATLAB/DFlow/influences_chap2/data/'...
    'fp_nolevees_elevation_1000m.mat']);

% Extract peak WSE from obs_fp_cell
max_km = max_dist/1000;
peak_WSE_1000m = zeros(max_km,6);
WSE_minus_bank = zeros(max_km,6);
WSE_minus_fp = zeros(max_km,6);
for i = 1:6
    a = 0;
    i_peak = obs_fp_cell{i,4};
    for j = 4:4:max_xs_num
        a = a + 1;
        peak_WSE_1000m(a,i) = mean(i_peak(j-3:j,3));
    end
    WSE_minus_bank(:,i) = peak_WSE_1000m(:,i) - ...
        bank_elevation_1000m(1:max_km,1);
    WSE_minus_fp(:,i) = peak_WSE_1000m(:,i) - ...
        fp_nolevees_elevation_1000m(1:max_km,1);
end

%% Plot QAR rates vs WSE_minus_bank and WSE_minus_fp

% Load QAR rates
load(['/Users/Colin/Documents/MATLAB/DFlow/influences_chap2/data/'...
    'dist_qar_rate1000m.mat']);

figure(1)
set(gcf,'color','white')
subplot(2,1,1)
plot(WSE_minus_bank(:,1),dist_qar_rate1000m(:,1),'o')
hold on
plot(WSE_minus_bank(:,2),dist_qar_rate1000m(:,2),'x')
lsline
ylim([0 0.04]);
xlabel('Peak WSE - BZ');
ylabel('Delta QAR');
mdl = fitlm(WSE_minus_bank(:,1),dist_qar_rate1000m(:,1));
r_squared_ordinary(1,1) = mdl.Rsquared.Ordinary;
r_squared_adjusted(1,1) = mdl.Rsquared.Adjusted;
mdl = fitlm(WSE_minus_bank(:,2),dist_qar_rate1000m(:,2));
r_squared_ordinary(1,2) = mdl.Rsquared.Ordinary;
r_squared_adjusted(1,2) = mdl.Rsquared.Adjusted;
set(gca,'Fontsize',14);
% ax = gca;
% ax.YScale = 'log';

subplot(2,1,2)
plot(abs(WSE_minus_fp(:,1)),dist_qar_rate1000m(:,1),'o')
hold on
plot(abs(WSE_minus_fp(:,2)),dist_qar_rate1000m(:,2),'x')
lsline
ylim([0 0.04]);
xlabel('Peak WSE - FPZ');
ylabel('Delta QAR');
mdl = fitlm(abs(WSE_minus_fp(:,1)),dist_qar_rate1000m(:,1));
r_squared_ordinary(2,1) = mdl.Rsquared.Ordinary;
r_squared_adjusted(2,1) = mdl.Rsquared.Adjusted;
mdl = fitlm(abs(WSE_minus_fp(:,2)),dist_qar_rate1000m(:,2));
r_squared_ordinary(2,2) = mdl.Rsquared.Ordinary;
r_squared_adjusted(2,2) = mdl.Rsquared.Adjusted;
set(gca,'Fontsize',14);
% ax = gca;
% ax.YScale = 'log';

figure(2)
set(gcf,'color','white')
subplot(2,1,1)
semilogy(WSE_minus_bank(:,1),dist_qar_rate1000m(:,1),'o')
hold on
semilogy(WSE_minus_bank(:,2),dist_qar_rate1000m(:,2),'x')
P1 = polyfit(WSE_minus_bank(:,1), log(dist_qar_rate1000m(:,1)), 1);
fit = exp(polyval(P1, WSE_minus_bank(:,1)));
plot(WSE_minus_bank(:,1),fit,'b-')
P2 = polyfit(WSE_minus_bank(:,2), log(dist_qar_rate1000m(:,2)), 1);
fit = exp(polyval(P2, WSE_minus_bank(:,2)));
plot(WSE_minus_bank(:,2),fit,'r-')
%lsline
ylim([0 0.04]);
xlabel('Peak WSE - BZ');
ylabel('Delta QAR');
% mdl = fitlm(WSE_minus_bank(:,1),dist_qar_rate1000m(:,1));
% r_squared_ordinary(1,1) = mdl.Rsquared.Ordinary;
% r_squared_adjusted(1,1) = mdl.Rsquared.Adjusted;
% mdl = fitlm(WSE_minus_bank(:,2),dist_qar_rate1000m(:,2));
% r_squared_ordinary(1,2) = mdl.Rsquared.Ordinary;
% r_squared_adjusted(1,2) = mdl.Rsquared.Adjusted;
set(gca,'Fontsize',14);
% ax = gca;
% ax.YScale = 'log';

subplot(2,1,2)
semilogy(abs(WSE_minus_fp(:,1)),dist_qar_rate1000m(:,1),'o')
hold on
semilogy(abs(WSE_minus_fp(:,2)),dist_qar_rate1000m(:,2),'x')
P3 = polyfit(abs(WSE_minus_fp(:,1)), log(dist_qar_rate1000m(:,1)), 1);
fit = exp(polyval(P3, abs(WSE_minus_fp(:,1))));
plot(abs(WSE_minus_fp(:,1)),fit,'b-')
P4 = polyfit(abs(WSE_minus_fp(:,2)), log(dist_qar_rate1000m(:,2)), 1);
fit = exp(polyval(P4, abs(WSE_minus_fp(:,2))));
plot(abs(WSE_minus_fp(:,2)),fit,'r-')
%lsline
ylim([0 0.04]);
xlabel('Peak WSE - FPZ');
ylabel('Delta QAR');
% mdl = fitlm(abs(WSE_minus_fp(:,1)),dist_qar_rate1000m(:,1));
% r_squared_ordinary(2,1) = mdl.Rsquared.Ordinary;
% r_squared_adjusted(2,1) = mdl.Rsquared.Adjusted;
% mdl = fitlm(abs(WSE_minus_fp(:,2)),dist_qar_rate1000m(:,2));
% r_squared_ordinary(2,2) = mdl.Rsquared.Ordinary;
% r_squared_adjusted(2,2) = mdl.Rsquared.Adjusted;
set(gca,'Fontsize',14);
% ax = gca;
% ax.YScale = 'log';