% Calculates skewness for upstream and 18-km cross-sections after
% his_extract_chap1 has been run

%%%%%%%%%%%%%%
%                sum(delta_Qneg_i / delta_t) / n
% skewness = 1 - -------------------------------
%                sum(delta_Qpos_i / delta_t) / n
%%%%%%%%%%%%%%
% removes lowest X% of flows

fraction_ignore = 0.01;

%% Extract 18-km flow data from each model condition and flow scenario

c2010_18kmflow = cell(6,1);
hist_18kmflow = cell(6,1); 
p2010_18kmflow = cell(6,1);
l2010_18kmflow = cell(6,1);

for i = 1:6
    c2010flow_i = xs_2010_cell{i,4};
    c2010_18kmflow(i,1) = {c2010flow_i(:,72)};
    
    histflow_i = xs_hist_cell{i,4};
    hist_18kmflow(i,1) = {histflow_i(:,72)};
    
    p2010flow_i = xs_prelidar_cell{i,4};
    p2010_18kmflow(i,1) = {p2010flow_i(:,72)};
    
    l2010flow_i = xs_2010lidar_cell{i,4};
    l2010_18kmflow(i,1) = {l2010flow_i(:,72)};
end

clear c2010flow_i histflow_i p2010flow_i l2010flow_i

%% Call function to extract skewness of upstream hydrographs

fpath_upstream = ...
    ['/Users/Colin/Documents/MATLAB/DFlow/roughness_investigation/'...
    'floods/'];
s180name = 's180peak.csv';
s300name = 's300peak.csv';
m180name = 'm180peak.csv';
m300name = 'm300peak.csv';
l180name = 'l180peak.csv';
l300name = 'l300peak.csv';

upstream(1,1) = {csvread(strcat(fpath_upstream,s180name))};
upstream(2,1) = {csvread(strcat(fpath_upstream,s300name))};
upstream(3,1) = {csvread(strcat(fpath_upstream,m180name))};
upstream(4,1) = {csvread(strcat(fpath_upstream,m300name))};
upstream(5,1) = {csvread(strcat(fpath_upstream,l180name))};
upstream(6,1) = {csvread(strcat(fpath_upstream,l300name))};

upstream_skew = zeros(6,1);

for i = 1:6
    
    up_flow = upstream{i,1};
    up_flow_seconds = up_flow(:,1)*60;
    interp_secs = 0:600:(5760*60);
    upstream_Q = interp1(up_flow_seconds,up_flow(:,2),interp_secs)';
    upstream_skew(i,1) = skewness_extract(upstream_Q,...
        max_Q(i,1)-150, fraction_ignore);
    
end

%% Call function to extract skewness of downstream hydrographs

hist_skew = zeros(6,1);
c2010_skew = zeros(6,1);
p2010_skew = zeros(6,1);
l2010_skew = zeros(6,1);

for j = 1:6
    
    hist_skew(j,1) = skewness_extract(hist_18kmflow{j,1},...
        max_Q(j,1)-150, fraction_ignore);
    c2010_skew(j,1) = skewness_extract(c2010_18kmflow{j,1},...
        max_Q(j,1)-150, fraction_ignore);
    p2010_skew(j,1) = skewness_extract(p2010_18kmflow{j,1},...
        max_Q(j,1)-150, fraction_ignore);
    l2010_skew(j,1) = skewness_extract(l2010_18kmflow{j,1},...
        max_Q(j,1)-150, fraction_ignore);
    
end

%% Call function to determine change in skewness from upstream

hist_delta_skew = zeros(6,1);
c2010_delta_skew = zeros(6,1);
p2010_delta_skew = zeros(6,1);
l2010_delta_skew = zeros(6,1);

for j = 1:6
    
    hist_delta_skew(j,1) = hist_skew(j,1) - upstream_skew(j,1);
    c2010_delta_skew(j,1) = c2010_skew(j,1) - upstream_skew(j,1);
    p2010_delta_skew(j,1) = p2010_skew(j,1) - upstream_skew(j,1);
    l2010_delta_skew(j,1) = l2010_skew(j,1) - upstream_skew(j,1);
    
end

