% Plots positive and negative fluxes for each interface type in comparison
% to total positive and negative flux values

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

%% Import model results of interest
% Currently set to use his_extract.m cell array output

[num_runs,~] = size(distance_cell);

for a = 1:num_runs
    
    run_row = a;
    
    lb_Ma_hold = Ma_cell{run_row,4};
    rb_Ma_hold = Ma_cell{run_row,5};
    lb_Mo_hold = Mo_cell{run_row,6};
    rb_Mo_hold = Mo_cell{run_row,7};
    
    [num_lb,~] = size(lb_Ma_hold);
    [num_rb,~] = size(rb_Ma_hold);
    
    if a == 1
        
        lb_xy = distance_cell{run_row,2};
        rb_xy = distance_cell{run_row,3};
        
        type_Ma_pos_neg = cell(2,num_runs);
        type_Mo_pos_neg = cell(2,num_runs);
        
        %% Import polygons associated with type of interface

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

        % Import inset floodplain polygons
        inset_interface = importdata(...
            ['/Users/Colin/Documents/MATLAB/DFlow/'...
            'mass_Mo_chap3/abq_inset_fp_pts']);

        %% Determine which type of flooplain model results apply to

        inset_num = max(inset_interface(:,3));
        rest_num = max(rest_interface(:,3));

        lb_type = zeros(num_lb,1);
        rb_type = zeros(num_rb,1);

        % Check if element is at inset floodplain interface
        for j = 1:inset_num+1
            inset_lb_j_yn = inpolygon(lb_xy(:,1), lb_xy(:,2),...
                inset_interface(inset_interface(:,3)==j-1,1),...
                inset_interface(inset_interface(:,3)==j-1,2));
            inset_rb_j_yn = inpolygon(rb_xy(:,1), rb_xy(:,2),...
                inset_interface(inset_interface(:,3)==j-1,1),...
                inset_interface(inset_interface(:,3)==j-1,2));
            lb_type(inset_lb_j_yn==1,1) = 1;
            rb_type(inset_rb_j_yn==1,1) = 1;
        end

        % Check if element is at restoration interface
        for k = 1:rest_num+1
            rest_lb_k_yn = inpolygon(lb_xy(:,1), lb_xy(:,2),...
                rest_interface(rest_interface(:,3)==k-1,1),...
                rest_interface(rest_interface(:,3)==k-1,2));
            rest_rb_k_yn = inpolygon(rb_xy(:,1), rb_xy(:,2),...
                rest_interface(rest_interface(:,3)==k-1,1),...
                rest_interface(rest_interface(:,3)==k-1,2));
            lb_type(rest_lb_k_yn==1,1) = 2;
            rb_type(rest_rb_k_yn==1,1) = 2;
        end
    
    end
    
    % Combine left and right interfaces
    all_types = [lb_type; rb_type];
    all_Ma_hold = [lb_Ma_hold; rb_Ma_hold];
    all_Mo_hold = [lb_Mo_hold; rb_Mo_hold];
    
    % Calculate positive and negative fluxes for all interfaces and each
    % type
    type_pos_Ma = zeros(4,t_steps);
    type_neg_Ma = zeros(4,t_steps);
    type_pos_Mo = zeros(4,t_steps);
    type_neg_Mo = zeros(4,t_steps);
    
    for n = 1:t_steps
        % Positive mass sums
        type_pos_Ma(4,n) = sum(all_Ma_hold(all_Ma_hold(:,n)>0,n));
        type_pos_Ma(1,n) = sum(all_Ma_hold(...
            all_Ma_hold(:,n)>0 & all_types==0,n));
        type_pos_Ma(2,n) = sum(all_Ma_hold(...
            all_Ma_hold(:,n)>0 & all_types==1,n));
        type_pos_Ma(3,n) = sum(all_Ma_hold(...
            all_Ma_hold(:,n)>0 & all_types==2,n));
        % Negative mass sums
        type_neg_Ma(4,n) = sum(all_Ma_hold(all_Ma_hold(:,n)<0,n));
        type_neg_Ma(1,n) = sum(all_Ma_hold(...
            all_Ma_hold(:,n)<0 & all_types==0,n));
        type_neg_Ma(2,n) = sum(all_Ma_hold(...
            all_Ma_hold(:,n)<0 & all_types==1,n));
        type_neg_Ma(3,n) = sum(all_Ma_hold(...
            all_Ma_hold(:,n)<0 & all_types==2,n));
        % Positive momentum sums
        type_pos_Mo(4,n) = sum(all_Mo_hold(all_Mo_hold(:,n)>0,n));
        type_pos_Mo(1,n) = sum(all_Mo_hold(...
            all_Mo_hold(:,n)>0 & all_types==0,n));
        type_pos_Mo(2,n) = sum(all_Mo_hold(...
            all_Mo_hold(:,n)>0 & all_types==1,n));
        type_pos_Mo(3,n) = sum(all_Mo_hold(...
            all_Mo_hold(:,n)>0 & all_types==2,n));
        % Negative momentum sums
        type_neg_Mo(4,n) = sum(all_Mo_hold(all_Mo_hold(:,n)<0,n));
        type_neg_Mo(1,n) = sum(all_Mo_hold(...
            all_Mo_hold(:,n)<0 & all_types==0,n));
        type_neg_Mo(2,n) = sum(all_Mo_hold(...
            all_Mo_hold(:,n)<0 & all_types==1,n));
        type_neg_Mo(3,n) = sum(all_Mo_hold(...
            all_Mo_hold(:,n)<0 & all_types==2,n));
    end
    
    type_Ma_pos_neg(1,a) = {type_pos_Ma};
    type_Ma_pos_neg(2,a) = {type_neg_Ma};
    
    type_Mo_pos_neg(1,a) = {type_pos_Mo};
    type_Mo_pos_neg(2,a) = {type_neg_Mo};
    
    [max_Ma_pos,max_Ma_pos_I] = max(type_pos_Ma(4,:));
    [max_Ma_neg,max_Ma_neg_I] = min(type_neg_Ma(4,:));
    [max_Mo_pos,max_Mo_pos_I] = max(type_pos_Mo(4,:));
    [max_Mo_neg,max_Mo_neg_I] = min(type_neg_Mo(4,:));
    
    time_max_Ma(a,1) = model_time(max_Ma_pos_I,1)/3600;
    time_max_Ma(a,2) = model_time(max_Ma_neg_I,1)/3600;
    time_max_Mo(a,1) = model_time(max_Mo_pos_I,1)/3600;
    time_max_Mo(a,2) = model_time(max_Mo_neg_I,1)/3600;
    
    
end

%% Plot mass and momentum flux by type for medium duration events

h1 = figure(1);
set(gcf,'Color','white');
model_time_hr = model_time / 3600;
for b = 3:4
    Ma_pos_plot = type_Ma_pos_neg{1,b};
    Ma_neg_plot = type_Ma_pos_neg{2,b};
    
    Mo_pos_plot = type_Mo_pos_neg{1,b};
    Mo_neg_plot = type_Mo_pos_neg{2,b};
    
    sp1 = subplot(2,2,b-2);
    plot(model_time_hr,Ma_pos_plot(4,:),'k-');
    hold on
    plot(model_time_hr,Ma_pos_plot(1,:));
    ax = gca;
    ax.ColorOrderIndex = 5;
    plot(model_time_hr,Ma_pos_plot(2,:));
    ax = gca;
    ax.ColorOrderIndex = 2;
    plot(model_time_hr,Ma_pos_plot(3,:));
    ax = gca;
    ax.ColorOrderIndex = 1;
    plot(model_time_hr,Ma_neg_plot(4,:),'k--');
    plot(model_time_hr,Ma_neg_plot(1,:),'--');
    ax = gca;
    ax.ColorOrderIndex = 5;
    plot(model_time_hr,Ma_neg_plot(2,:),'--');
    ax = gca;
    ax.ColorOrderIndex = 2;
    plot(model_time_hr,Ma_neg_plot(3,:),'--');
    xlim([0 40])
    %ylim([0 .075])
    %set(gca,'Ytick',[0 0.025 0.05 0.075]);
    if b == 3;
         ylabel('Mass Flux (kg/s)');
         %xlabel('Probability');
    elseif b == 4;
         %xlabel('Time (hr)');
         %set(gca,'yticklabel','');
    end
    set(gca,'Fontsize',14);
    box(sp1,'off')
    hold off

    sp2=subplot(2,2,b);
    plot(model_time_hr,Mo_pos_plot(4,:),'k-');
    hold on
    plot(model_time_hr,Mo_pos_plot(1,:));
    ax = gca;
    ax.ColorOrderIndex = 5;
    plot(model_time_hr,Mo_pos_plot(2,:));
    ax = gca;
    ax.ColorOrderIndex = 2;
    plot(model_time_hr,Mo_pos_plot(3,:));
    ax = gca;
    ax.ColorOrderIndex = 1;
    plot(model_time_hr,Mo_neg_plot(4,:),'k--');
    ax = gca;
    ax.ColorOrderIndex = 1;
    plot(model_time_hr,Mo_neg_plot(1,:),'--');
    ax = gca;
    ax.ColorOrderIndex = 5;
    plot(model_time_hr,Mo_neg_plot(2,:),'--');
    ax = gca;
    ax.ColorOrderIndex = 2;
    plot(model_time_hr,Mo_neg_plot(3,:),'--');
    xlim([0 40])
    %ylim([0 .075])
    %set(gca,'Ytick',[0 0.025 0.05 0.075]);
    if b == 3;
         %set(gca,'xticklabel','');
         xlabel('Time (hr)');
         ylabel('Momentum Flux (N)');
    elseif b == 4;
         %set(gca,'xticklabel','');
         %set(gca,'yticklabel','');
         %ylabel('Momentum Flux (N)');
         xlabel('Time (hr)');
    end
    set(gca,'Fontsize',14);
    box(sp2,'off')
    hold off
end

