% Computes attenuation of distributed runs

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'));

% Collect baseflow maximum for distributed FP n

base30cmsdist_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3051_30cmsdist/'...
    'abq10a12bf30cmsdistrst01_0000_20160104_090000_his.nc'];
base150cmsdist_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3052_150cmsdist/'...
    'abq10a12bf150cmsdistrst01_0000_20160102_230000_his.nc'];

% Extract flow at all cross-section times
base30cmsdist_all_Q = ...
    ncread(base30cmsdist_file,'cross_section_discharge');
base150cmsdist_all_Q = ...
    ncread(base150cmsdist_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
base30cmsdist_Q = base30cmsdist_all_Q(:,end);
base150cmsdist_Q = base150cmsdist_all_Q(:,end);

clear base30cmsdist_all_Q base150cmsdist_all_Q

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

peak180cmsdist_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3083_180cmsbank040/'...
    'abq10a12bf180cmsbank040rst02_0000_20160104_160000_his.nc'];
peak300cmsdist_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3046_300cmsdist/'...
    'abq10a12bf300cmsrst02_0000_20160103_210000_his.nc'];

% Extract flow at all cross-section times
peak180cmsdist_all_Q = ncread(peak180cmsdist_file,'cross_section_discharge');
peak300cmsdist_all_Q = ncread(peak300cmsdist_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
peak180cmsdist_Q = peak180cmsdist_all_Q(:,end);
peak300cmsdist_Q = peak300cmsdist_all_Q(:,end);

clear peak180cmsdist_all_Q peak300cmsdist_all_Q

% Collect peak flow information from distributed FP n 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));
xs_fpdist_cell = cell(folder_num,7);

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

    if mod(fi,2) == 0
        xs_peak_Q = peak300cmsdist_Q;
        xs_base_Q = base150cmsdist_Q;
    elseif mod(fi,2) == 1
        xs_peak_Q = peak180cmsdist_Q;
        xs_base_Q = base30cmsdist_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_fpdist_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

%% Calculate rates of attenuation

max_xs_id = max_dist / 250;
dist_qar_rates = zeros(max_xs_id-1,6);
dist_qar_proportional = zeros(max_xs_id-1,6);

for i = 1:6
    
    flood_qar = xs_fpdist_cell{i,9};
    total_qar = 1 - flood_qar(max_xs_id,1);
    
    for j = 2:max_xs_id
        
        dist_qar_rates(j-1,i) = flood_qar(j-1,1) - flood_qar(j,1);
        dist_qar_proportional(j-1,i) = dist_qar_rates(j-1,i) / total_qar;
        
    end
    
end