% Analyzes discharge fluxes for kilometer cross-sections and
% channel-floodplain interfaces to understand lateral and longitudinal flow
% relationships

clearvars -except chan_dist_max distance_cell Ma_cell max_min_cell ...
    max_t max_t_r min_t min_t_r Mo_cell model_time Q_cell raw_cell ...
    segment250pts t_steps xs_cell

%% Determine which interface elements belong to which kilometer

max_xs = chan_dist_max / 250;
km_max = chan_dist_max/1000;
[num_runs,~] = size(distance_cell);

lb_xy = distance_cell{1,2};
rb_xy = distance_cell{1,3};
all_xy = [lb_xy; rb_xy];

km_id = zeros((numel(lb_xy(:,1))+numel(rb_xy(:,1))),1);
for i = 4:4:max_xs
    % Run inpolygon for previous three 250 m segments
    for j = 1:3
        seg250_yn = inpolygon(all_xy(:,1),all_xy(:,2),...
            segment250pts(segment250pts(:,3)==(i-j),1),...
            segment250pts(segment250pts(:,3)==(i-j),2));
        km_id(seg250_yn==1,1) = i;
    end
    % Run inpolygon for current 250 m segments
    seg250_yn = inpolygon(all_xy(:,1),all_xy(:,2),...
    segment250pts(segment250pts(:,3)==(i),1),...
    segment250pts(segment250pts(:,3)==(i),2));
    km_id(seg250_yn==1,1) = i/4;
end

%% Calculate lateral, longitudinal exchange and ratios

lateral_Q = cell(6,1);
longitudinal_Q = cell(6,1);
lat_long_ratios = cell(6,1);

for k = 1:num_runs
    
    % Calculate lateral discharge within the km upstream of cross-section
    all_k_lat = [Q_cell{k,2}; Q_cell{k,3}];
    lat_segs_Q = zeros(km_max,t_steps);
    for l = 1:km_max
        lat_segs_Q(l,:) = sum(all_k_lat(km_id==l,:));
    end
    lateral_Q(k,1) = {lat_segs_Q};

    % Calculate longitudinal discharge at each cross-section
    all_k_long = [xs_cell{k,4}];
    xs_dists = [xs_cell{k,6}];
    long_segs_Q = zeros(km_max,t_steps);
    for m = 1:km_max
        long_segs_Q(m,:) = all_k_long(:,xs_dists(:,1)==m*1000)';
    end
    longitudinal_Q(k,1) = {long_segs_Q};

    % Calculate lateral-longitudinal ratios
    ratios_n = zeros(km_max,t_steps);
    for n = 1:km_max
        ratios_n(n,:) = lat_segs_Q(n,:) ./ long_segs_Q(n,:);
    end
    lat_long_ratios(k,1) = {ratios_n};
    
    ratios_n_mean(:,k) = mean(ratios_n,2);
    ratios_n_max(:,k) = max(ratios_n,[],2);
    ratios_n_min(:,k) = min(ratios_n,[],2);
    
    ratios_n_max_percent(:,k) = sum(ratios_n_max(:,k)<0.1) / ...
        numel(ratios_n_max(:,k));
    
    clear ratios_n
end

ratios_less_tenth = mean(ratios_n_max_percent);

%% Plot ratios

% Plot m180 flood peak for starters
%ratio_m180p = lat_long_ratios{3,1};

f1=figure(1);
set(f1,'Color','white');
fsize=14;
spwidth = 0.40;
spheight = 0.25;
for o = 1:num_runs
    ratio_o = lat_long_ratios{o,1};
    hand1 = subplot(3,2,o);
    phand = get(hand1,'pos');
    for p = 11:15%4:km_max
        plot(model_time(min_t_r:max_t_r,1)/3600,...
            ratio_o(p,min_t_r:max_t_r));
        hold on
    end
    if o == 1
        legend('11 km','12 km','13 km','14 km','15 km');
        legend('Location','east')
        set(gca,'xticklabel','');
        phand(2) = 0.71;
        phand(3) = spwidth;
        phand(4) = spheight;
        %title('Low baseflow')
    elseif o == 2;
        set(gca,'yticklabel','');
        set(gca,'xticklabel','');
        phand(2) = 0.71;
        phand(3) = spwidth;
        phand(4) = spheight;
        %title('High baseflow')
    elseif o == 3
        set(gca,'xticklabel','');
        ylabel('Ratio of lateral to longitudinal discharge')
        phand(2) = 0.41;
        phand(3) = spwidth;
        phand(4) = spheight;
    elseif o == 4;
        set(gca,'yticklabel','');
        set(gca,'xticklabel','');
        phand(2) = 0.41;
        phand(3) = spwidth;
        phand(4) = spheight;
    elseif o == 5;
        xlabel('Time (hr)');
        phand(2) = 0.11;
        phand(3) = spwidth;
        phand(4) = spheight;
    elseif o == 6;
        set(gca,'yticklabel','');
        xlabel('Time (hr)');
        phand(2) = 0.11;
        phand(3) = spwidth;
        phand(4) = spheight;
    end
    set(hand1, 'pos', phand);
    xlim([0 80])
    ylim([-0.1 0.2])
    set(gca,'Fontsize',fsize);
    hold off
end
set(findall(f1, 'Type', 'Text'),'FontWeight', 'Normal')