% Extract information from *_his.nc files

clear

%% INPUTS

chan_dist_max = 27000; % km
max_Q = [180; 300; 180; 300; 180; 300]; % m^3/s
site_num = 15; % number of sites
min_t = 7200; % plots start at 2 hour
max_t = 96 * 3600; % hours * seconds in an hour - gives seconds
fpath_root = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/'];

% Define folder with *his.nc files
fpath = ['run3034_s180pback/';'run3035_s300pback/';...
    'run3036_m180pback/';'run3037_m300pback/';...
    'run3038_l180pback/';'run3039_l300pback/']; % file path
% files = struct2cell(dir(fpath)); % Get files from folder
% [~,file_count] = size(files);
% fnames = files(1,3:file_count)'; % extract file names
% [fnum,~] = size(fnames); % number of files

segment250pts = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'mass_Mo_chap3/abq_250m_segment_polygons'],' '); % import 250 segments

folder_num = numel(fpath(:,1));

xs_cell = cell(folder_num,7);
raw_cell = cell(folder_num,3);
distance_cell = cell(folder_num,9);
Q_cell = cell(folder_num,7);
Ma_cell = cell(folder_num,13);
Mo_cell = cell(folder_num,15);
max_min_cell = cell(folder_num,37);

for fi = 1:folder_num
    
    xs_cell(fi,1) = {fpath(fi,:)};
    raw_cell(fi,1) = {fpath(fi,:)};
    distance_cell(fi,1) = {fpath(fi,:)};
    Q_cell(fi,1) = {fpath(fi,:)};
    Ma_cell(fi,1) = {fpath(fi,:)};
    Mo_cell(fi,1) = {fpath(fi,:)};
    max_min_cell(fi,1) = {fpath(fi,:)};
    
    max_Q_fi = max_Q(fi,1);
    
    fpath_fi = strcat(fpath_root,fpath(fi,:));
    
    fnames = struct2cell(dir(strcat(fpath_fi,'*his.nc')))'; 
    [fnum,~] = size(fnames); % number of files

    % Import channel/island areas
    load(['/Users/Colin/Documents/MATLAB/DFlow/mass_Mo_chap3/'...
        'rg_channel_island_area.mat']);

    % Import restoration site connection polygons
    rest_con = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
        'mass_Mo_chap3/abq12_restoration_feature_connections_pts']);

    rs_index = cell(1,site_num);
    % Observation ID array
    for s = 1:site_num
        rs_index(s) = {strcat('s',num2str(s))}; % create cell with site ids
    end
    obs_marker = ['l', 'r']; % array marking left or right boundaries

    %% 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
    
    [xs_cell(fi,2:7), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his(fnum, fnames, fpath_fi, chan_dist_max, min_t, ...
        max_t, max_Q_fi);

    %% Extract channel-floodplain boundary information (separate elements)
    % raw_cell
    % raw_cell columns (label, left, right)
    % Cells within left and right cells
    % Row 1 - element coordinates (NAD83 m)
    % Row 2 - Vx (m/s)
    % Row 3 - Vy (m/s)
    % Row 4 - Depth (m)
    % Row 5 - midpoint of flux line coordinate
    % Row 6 - Discharge (m^3/s)
    % Row 7 - Momentum flux per area (N/m^2)
    % Row 8 - Momentum flux (N)
    % Row 9 - Momentum length (N/m)
    % Row 10 - Upstream and downstream distances (m)
    % Row 11 - Cumulative downstream distance (m)
    % Row 12 - Area of flux (m^2)
    % Row 13 - Interface element length (m)

    [raw_cell(fi,2:3), model_time] = ...
        obs2fpboundaries(fnum,fpath_fi,fnames,obs_marker);
    
    %% Extract location and distance information in the following form
    % distance_cell
    % Column 1 - Run Label
    % Column 2 - left midpoint coordinates
    % Column 3 - right midpoint coordinates
    % Column 4 - left cumulative downstream distance
    % Column 5 - right cumulative downstream distance
    % Column 6 - left upstream and downstream distances from element
    % Column 7 - right upstream and downstream distances from element
    % Column 8 - left center-normalized distance downstream
    % Column 9 - right center-normalized distance downstream
    
    distance_cell(fi,2:9) = distance_calcs(raw_cell{fi,2},...
        raw_cell{fi,3}, segment250pts);

    %% Extract discharge data in the following form
    % Q_cell
    % Column 1 - Run Label
    % Column 2 - left raw discharge
    % Column 3 - right raw discharge
    % Column 4 - left baseflow discharge
    % Column 5 - right baseflow discharge
    % Column 6 - left above baseflow discharge
    % Column 7 - right above baseflow discharge
    % Column 8 - left temporal discharge along 27 km (above base)
    % Column 9 - right temporal discharge along 27 km (above base)
    
    Q_cell(fi,2:9) = Q_calcs(raw_cell{fi,2}, raw_cell{fi,3}, t_steps);

    %% Extract mass flux information in the following form
    % Ma_cell
    % Column 1 - Run Label
    % Column 2 - left mass flux with baseflow
    % Column 3 - right mass flux with baseflow
    % Column 4 - left mass flux WITHOUT baseflow
    % Column 5 - right mass flux WITHOUT baseflow
    % Column 6 - left mass flux per unit length WITHOUT baseflow
    % Column 7 - right mass flux per unit length WITHOUT baseflow
    % Column 8 - left temporal mass flux along 27 km (above base)
    % Column 9 - right temporal mass flux along 27 km (above base)
    % Column 10 - left positive temporal flux along 27 km (above base)
    % Column 11 - left negative temporal flux along 27 km (above base)
    % Column 12 - right positive temporal flux along 27 km (above base)
    % Column 13 - right negative temporal flux along 27 km (above base)
    
    Ma_cell(fi,2:13) = Ma_calcs(Q_cell{fi,6}, Q_cell{fi,7},...
        Q_cell{fi,2}, Q_cell{fi,3}, t_steps, raw_cell{fi,2},...
        raw_cell{fi,3});

    %% Extract momentum flux information in the following form
    % Mo_cell
    % Column 1 - Run Label
    % Column 2 - left momentum flux with baseflow
    % Column 3 - right momentum flux with baseflow
    % Column 4 - left momentum flux AT baseflow
    % Column 5 - right momentum flux AT baseflow
    % Column 6 - left momentum flux WITHOUT baseflow
    % Column 7 - right momentum flux WITHOUT baseflow
    % Column 8 - left momentum flux per unit length WITHOUT baseflow
    % Column 9 - right momentum flux per unit length WITHOUT baseflow
    % Column 10 - left temporal momentum flux along 27 km (above base)
    % Column 11 - right temporal momentum flux along 27 km (above base)
    % Column 12 - left positive temporal flux along 27 km (above base)
    % Column 13 - left negative temporal flux along 27 km (above base)
    % Column 14 - right positive temporal flux along 27 km (above base)
    % Column 15 - right negative temporal flux along 27 km (above base)
    
    Mo_cell(fi,2:15) = Mo_calcs(raw_cell{fi,2}, raw_cell{fi,3}, t_steps);

    %% Maximum and minimum discharge, mass, and momentum extraction
    % Mo_cell
    % Column 1 - Run Label
    % Column 2 - left Q max and time index above baseflow
    % Column 3 - right Q max and time index above baseflow
    % Column 4 - left Q min and time index above baseflow
    % Column 5 - right Q min and time index above baseflow
    % Column 6 - left Q max, index, location above baseflow (> 0.01 m^3/s)
    % Column 7 - right Q max, index, location above baseflow (> 0.01 m^3/s)
    % Column 8 - left Q min, index, location above baseflow (> 0.01 m^3/s)
    % Column 9 - right Q min, index, location above baseflow (> 0.01 m^3/s)
    % Column 10 - left Q max, index, location above baseflow (> 0.01 m^3/s)
    %               - sorted
    % Column 11 - right Q max, index, location above baseflow (> 0.01
    %               m^3/s) - sorted
    % Column 12 - left Q min, index, location above baseflow (> 0.01 m^3/s)
    %               - sorted
    % Column 13 - right Q min, index, location above baseflow (> 0.01
    %               m^3/s) - sorted
    % Columns 14 - 25 - Reapeat Column 2-13 for mass not Q
    % Columns 26 - 37 - Reapeat Column 2-13 for momentum not Q
    
    [max_min_cell(fi,2:37)] = max_min_calcs(distance_cell(fi,:),...
        Q_cell(fi,:), Ma_cell(fi,:), Mo_cell(fi,:));

end

%% Clear variables
clear fi fnames fnum folder_num fpath fpath_fi fpath_root max_Q ...
    max_Q_fi obs_marker rest_con rg_channel_island_area rs_index s site_num

%% Split mass and momentum into 250m segments
% xs_segs = zeros(max_xs_num,5); 
% 
% for xs_i = 1:max_xs_num
%     xs_segs(xs_i,1) = xs_i;
%     segment_i = segment250pts(segment250pts(:,3)==xs_i,:);
%     timeofmax_i = xs_max_Q_27km(xs_i,2);
%     lin_poly_yn = inpolygon(lb_xy(:,1),lb_xy(:,2),...
%         segment_i(:,1),segment_i(:,2));
%     rin_poly_yn = inpolygon(rb_xy(:,1),rb_xy(:,2),...
%         segment_i(:,1),segment_i(:,2));
%     xs_segs(xs_i,2) = sum(lin_poly_yn(:)) + sum(rin_poly_yn(:));
%     xs_segs(xs_i,3) = sum(lb_Ma(lin_poly_yn(:)==1,timeofmax_i/600)) +...
%         sum(rb_Ma(rin_poly_yn(:)==1,timeofmax_i/600));
%     xs_segs(xs_i,4) = sum(lb_Moflux(lin_poly_yn(:)==1,timeofmax_i/600)) ...
%         + sum(rb_Moflux(rin_poly_yn(:)==1,timeofmax_i/600));
%     clear segment_i
% end
% 
% for xs_i = 1:max_xs_num
%     if xs_i == 1
%         xs_segs(xs_i,5) = xs_segs(xs_i,3);
%         xs_segs(xs_i,6) = xs_segs(xs_i,4);
%     else
%         xs_segs(xs_i,5) = xs_segs(xs_i,3) + xs_segs(xs_i-1,5);
%         xs_segs(xs_i,6) = xs_segs(xs_i,4) + xs_segs(xs_i-1,6);
%     end
% end
% 
% %% Restoration connection extraction
% [rest_con_r,~] = size(rest_con);
% rc_data = cell(7,3);
% 
% for rest_i = 10:10:70
%     lupcon_yn = inpolygon(lb_xy(:,1),lb_xy(:,2),...
%         rest_con(rest_con(:,3)==rest_i,1),...
%         rest_con(rest_con(:,3)==rest_i,2));
%     ldowncon_yn = inpolygon(lb_xy(:,1),lb_xy(:,2),...
%         rest_con(rest_con(:,3)==rest_i+1,1),...
%         rest_con(rest_con(:,3)==rest_i+1,2));
%     rupcon_yn = inpolygon(rb_xy(:,1),rb_xy(:,2),...
%         rest_con(rest_con(:,3)==rest_i,1),...
%         rest_con(rest_con(:,3)==rest_i,2));
%     rdowncon_yn = inpolygon(rb_xy(:,1),rb_xy(:,2),...
%         rest_con(rest_con(:,3)==rest_i+1,1),...
%         rest_con(rest_con(:,3)==rest_i+1,2));
%     
%     lrestcon_yn = lupcon_yn + ldowncon_yn;
%     rrestcon_yn = rupcon_yn + rdowncon_yn;
%     
%     rc_data{rest_i/10,1} = rest_i/10;
%     
%     if sum(lrestcon_yn) > 0
%         rc_data{rest_i/10,2} = lb_Q(lrestcon_yn>0,:);
%         rc_data{rest_i/10,3} = lb_Ma(lrestcon_yn>0,:);
%         rc_Ma = lb_Ma(lrestcon_yn>0,:);
%         rc_Ma_pos = zeros(1,t_steps);
%         rc_Ma_neg = zeros(1,t_steps);
%         for rc_i = 1:t_steps
%             rc_Ma_pos(1,rc_i) =...
%                 sum(rc_Ma(rc_Ma(:,rc_i)>0,rc_i));
%             rc_Ma_neg(1,rc_i) =...
%                 sum(rc_Ma(rc_Ma(:,rc_i)<0,rc_i));
%         end
%         rc_data{rest_i/10,4} = rc_Ma_pos;
%         rc_data{rest_i/10,5} = rc_Ma_neg;
%     elseif sum(rrestcon_yn) > 0
%         rc_data{rest_i/10,2} = rb_Q(rrestcon_yn>0,:);
%         rc_data{rest_i/10,3} = rb_Ma(rrestcon_yn>0,:);
%         rc_Ma = lb_Ma(lrestcon_yn>0,:);
%         rc_Ma_pos = zeros(1,t_steps);
%         rc_Ma_neg = zeros(1,t_steps);
%         for rc_i = 1:t_steps
%             rc_Ma_pos(1,rc_i) =...
%                 sum(rc_Ma(rc_Ma(:,rc_i)>0,rc_i));
%             rc_Ma_neg(1,rc_i) =...
%                 sum(rc_Ma(rc_Ma(:,rc_i)<0,rc_i));
%         end
%         rc_data{rest_i/10,4} = rc_Ma_pos;
%         rc_data{rest_i/10,5} = rc_Ma_neg;
%     end
% end
