% Loads mass transfer files

%clear

%% Load all *.mat files

fdir = '/Users/Colin/Documents/MATLAB/DFlow/mass_Mo_chap3/bndry_Ma_vs_t';
s = what(fdir); %change dir for your directory name 
matfiles = s.mat;
for a=1:numel(matfiles)
    load(char(matfiles(a)))
end

%% Create min and max time indices

min_t = 7200; % plots start at 2 hour
max_t = 96 * 3600; % hours * seconds in an hour - gives seconds
[min_t_r,~] = find(xs_time(:,1)==min_t); % index of set minimum time
[max_t_r,~] = find(xs_time(:,1)==max_t); % index of set maximum time

%% Plot all files
fsize = 14;
h = figure(1);
h.Units = 'inches';
%h.Position = [1 1 11 15];
h.PaperSize = [8.5 9];
set(gcf,'Color','white')
max_x = 84;
min_x = 0;
max_y = 60000;
min_y = -40000;

sp1 = subplot(4,2,1);
sp1f1=fill([9 xs_max_Q_27km_s180pback(108,2)/3600 ...
    xs_max_Q_27km_s180pback(108,2)/3600 9],...
    [-10000 -10000 40000 40000],[.9 .9 .9],'EdgeColor','none');
hold on
sp1p1 = plot([0 360000],[0 0],'k-','LineWidth',0.25);
plot(xs_time(min_t_r:max_t_r,1)/3600,...
    rb_Ma_t_s180pback(1,min_t_r:max_t_r)/10000,'r-','LineWidth',1.5);
plot(xs_time(min_t_r:max_t_r,1)/3600,...
    lb_Ma_t_s180pback(1,min_t_r:max_t_r)/10000,'b-','LineWidth',1.5)
% plot([xs_max_Q_27km_s180pback(108,2)/3600 ...
%     xs_max_Q_27km_s180pback(108,2)/3600],...
%     [-40000 80000],'k:','LineWidth',1.5)
% plot(xs_time(min_t_r:max_t_r,1)/3600,...
%     rb_Ma_t_pres180pback(1,min_t_r:max_t_r),'r--')
% plot(xs_time(min_t_r:max_t_r,1)/3600,...
%     lb_Ma_t_pres180pback(1,min_t_r:max_t_r),'b--')
% plot([xs_max_Q_27km_pres180pback(108,2)/3600 ...
%     xs_max_Q_27km_pres180pback(108,2)/3600],...
%     [-40000 80000],'k--','LineWidth',1.5)
set(gca,'Fontsize',fsize)
%legend('Right CH-FP','Left CH-FP',sprintf('Flood peak\nleaves reach'),...
%    'Location','Northeast')
%xlabel('Time (hrs)')
%ylabel({'Mass flux'; 'across interface (kg/s)'})
% axis([min_x max_x min_y max_y]);
% set(gca,'Ytick',[-30000 0 30000 60000]);
axis([min_x max_x min_y/10000 max_y/10000]);
set(gca,'Ytick',[-3 0 3 6]);
set(gca,'xticklabel',[])
ylabel({'Net';'mass flux';'(kg/s*10^4)'})
set(gca,'Layer','top')
%grid on
hold off

sp3 = subplot(4,2,3);
fill([16 xs_max_Q_27km_m180pback(108,2)/3600 ...
    xs_max_Q_27km_m180pback(108,2)/3600 16],...
    [-10000 -10000 40000 40000],[.9 .9 .9],'EdgeColor','none')
hold on
plot([0 360000],[0 0],'k-','LineWidth',0.25);
plot(xs_time(min_t_r:max_t_r,1)/3600,...
    rb_Ma_t_m180pback(1,min_t_r:max_t_r)/10000,'r-','LineWidth',1.5)
plot(xs_time(min_t_r:max_t_r,1)/3600,...
    lb_Ma_t_m180pback(1,min_t_r:max_t_r)/10000,'b-','LineWidth',1.5)
% plot([xs_max_Q_27km_m180pback(108,2)/3600 ...
%     xs_max_Q_27km_m180pback(108,2)/3600],...
%     [-40000 80000],'k:','LineWidth',1.5)
% Plot pre restoration conditions
% plot(xs_time(min_t_r:max_t_r,1)/3600,...
%     rb_Ma_t_prem180pback(1,min_t_r:max_t_r),'r--')
% plot(xs_time(min_t_r:max_t_r,1)/3600,...
%     lb_Ma_t_prem180pback(1,min_t_r:max_t_r),'b--')
% plot([xs_max_Q_27km_prem180pback(108,2)/3600 ...
%     xs_max_Q_27km_prem180pback(108,2)/3600],...
%     [-40000 80000],'k--','LineWidth',1.5)
set(gca,'Fontsize',fsize)
% legend('Right CH-FP','Left CH-FP','Flood peak leaves reach',...
%     'Location','Northeast')
%xlabel('Time (hrs)')
ylabel({'Net';'mass flux';'(kg/s*10^4)'})
% axis([min_x max_x min_y max_y]);
% set(gca,'Ytick',[-30000 0 30000 60000]);
axis([min_x max_x min_y/10000 max_y/10000]);
set(gca,'Ytick',[-3 0 3 6]);
set(gca,'xticklabel',[])
set(gca,'Layer','top')
%grid on
hold off

sp7 = subplot(4,2,7);
fill([30 xs_max_Q_27km_l180pback(108,2)/3600 ...
    xs_max_Q_27km_l180pback(108,2)/3600 30],...
    [-10000 -10000 40000 40000],[.9 .9 .9],'EdgeColor','none')
hold on
plot([0 360000],[0 0],'k-','LineWidth',0.25);
plot(xs_time(min_t_r:max_t_r,1)/3600,...
    rb_Ma_t_l180pback(1,min_t_r:max_t_r)/10000,'r-','LineWidth',1.5)
plot(xs_time(min_t_r:max_t_r,1)/3600,...
    lb_Ma_t_l180pback(1,min_t_r:max_t_r)/10000,'b-','LineWidth',1.5)
% plot([xs_max_Q_27km_l180pback(108,2)/3600 ...
%     xs_max_Q_27km_l180pback(108,2)/3600],...
%     [-40000 80000],'k:','LineWidth',1.5)
set(gca,'Fontsize',fsize)
% legend('Right CH-FP','Left CH-FP','Flood peak leaves reach',...
%     'Location','Northeast')
xlabel('Time (hrs)')
%ylabel({'Mass flux'; 'across interface (kg/s)'})
% axis([min_x max_x min_y max_y]);
% set(gca,'Ytick',[-30000 0 30000 60000]);
axis([min_x max_x min_y/10000 max_y/10000]);
set(gca,'Ytick',[-3 0 3 6]);
ylabel({'Net';'mass flux';'(kg/s*10^4)'})
set(gca,'Layer','top')
%grid on
hold off

sp2 = subplot(4,2,2);
fill([9 xs_max_Q_27km_s300pback(108,2)/3600 ...
    xs_max_Q_27km_s300pback(108,2)/3600 9],...
    [-10000 -10000 40000 40000],[.9 .9 .9],'EdgeColor','none')
hold on
plot([0 360000],[0 0],'k-','LineWidth',0.25);
plot(xs_time(min_t_r:max_t_r,1)/3600,...
    rb_Ma_t_s300pback(1,min_t_r:max_t_r)/10000,'r-','LineWidth',1.5)
plot(xs_time(min_t_r:max_t_r,1)/3600,...
    lb_Ma_t_s300pback(1,min_t_r:max_t_r)/10000,'b-','LineWidth',1.5)
% plot([xs_max_Q_27km_s300pback(108,2)/3600 ...
%     xs_max_Q_27km_s300pback(108,2)/3600],...
%     [-40000 80000],'k:','LineWidth',1.5)
% plot(xs_time(min_t_r:max_t_r,1)/3600,...
%     rb_Ma_t_pres300pback(1,min_t_r:max_t_r),'r--')
% hold on
% plot(xs_time(min_t_r:max_t_r,1)/3600,...
%     lb_Ma_t_pres300pback(1,min_t_r:max_t_r),'b--')
% plot([xs_max_Q_27km_pres300pback(108,2)/3600 ...
%     xs_max_Q_27km_pres300pback(108,2)/3600],...
%     [-40000 80000],'k:','LineWidth',1.5)
set(gca,'Fontsize',fsize)
% legend('Right CH-FP','Left CH-FP','Flood peak leaves reach',...
%     'Location','Northeast')
%xlabel('Time (hrs)')
%ylabel('Mass flux'; 'across interface (kg/s)')
% axis([min_x max_x min_y max_y]);
% set(gca,'Ytick',[-30000 0 30000 60000]);
axis([min_x max_x min_y/10000 max_y/10000]);
set(gca,'Ytick',[-3 0 3 6]);
set(gca,'xticklabel',[])
set(gca,'Layer','top')
%grid on
hold off

sp4 = subplot(4,2,4);
fill([16 xs_max_Q_27km_m300pback(108,2)/3600 ...
    xs_max_Q_27km_m300pback(108,2)/3600 16],...
    [-10000 -10000 40000 40000],[.9 .9 .9],'EdgeColor','none')
hold on
plot([0 360000],[0 0],'k-','LineWidth',0.25);
plot(xs_time(min_t_r:max_t_r,1)/3600,...
    rb_Ma_t_m300pback(1,min_t_r:max_t_r)/10000,'r-','LineWidth',1.5)
plot(xs_time(min_t_r:max_t_r,1)/3600,...
    lb_Ma_t_m300pback(1,min_t_r:max_t_r)/10000,'b-','LineWidth',1.5)
% plot([xs_max_Q_27km_m300pback(108,2)/3600 ...
%     xs_max_Q_27km_m300pback(108,2)/3600],...
%     [-40000 80000],'k:','LineWidth',1.5)
set(gca,'Fontsize',fsize)
% legend('Right CH-FP','Left CH-FP','Flood peak leaves reach',...
%     'Location','Northeast')
%xlabel('Time (hrs)')
%ylabel({'Mass flux'; 'across interface (kg/s)'})
% axis([min_x max_x min_y max_y]);
% set(gca,'Ytick',[-30000 0 30000 60000]);
axis([min_x max_x min_y/10000 max_y/10000]);
set(gca,'Ytick',[-3 0 3 6]);
set(gca,'xticklabel',[])
set(gca,'Layer','top')
%grid on
hold off

sp8 = subplot(4,2,8);
fill([30 xs_max_Q_27km_l300pback(108,2)/3600 ...
    xs_max_Q_27km_l300pback(108,2)/3600 30],...
    [-10000 -10000 40000 40000],[.9 .9 .9],'EdgeColor','none')
hold on
plot([0 360000],[0 0],'k-','LineWidth',0.25);
plot(xs_time(min_t_r:max_t_r,1)/3600,...
    rb_Ma_t_l300pback(1,min_t_r:max_t_r)/10000,'r-','LineWidth',1.5)
plot(xs_time(min_t_r:max_t_r,1)/3600,...
    lb_Ma_t_l300pback(1,min_t_r:max_t_r)/10000,'b-','LineWidth',1.5)
% plot([xs_max_Q_27km_l300pback(108,2)/3600 ...
%     xs_max_Q_27km_l300pback(108,2)/3600],...
%     [-40000 80000],'k:','LineWidth',1.5)
set(gca,'Fontsize',fsize)
% legend('Right CH-FP','Left CH-FP','Flood peak leaves reach',...
%     'Location','Northeast')
xlabel('Time (hrs)')
%ylabel('Mass flux'; 'across interface (kg/s)')
axis([min_x max_x min_y/10000 max_y/10000]);
set(gca,'Ytick',[-3 0 3 6]);
set(gca,'Layer','top')
%grid on
hold off

%% Add subplots for total fluxes

% 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 min_x max_x sp1 sp2 sp3 sp4 sp7 sp8 ...
    sp1p1 xs_max_Q_27km_m180pback xs_max_Q_27km_m300pback

%% 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};
        
end

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};
    
    if b == 3
        sp5 = subplot(4,2,b+2);
    elseif b == 4;
        sp6 = subplot(4,2,b+2);
    end
    if b == 3
        pfi1=fill([16 xs_max_Q_27km_m180pback(108,2)/3600 ...
            xs_max_Q_27km_m180pback(108,2)/3600 16],...
            [-10000 -10000 40000 40000],[.9 .9 .9],'EdgeColor','none');
    elseif b == 4
        pfi2=fill([16 xs_max_Q_27km_m300pback(108,2)/3600 ...
            xs_max_Q_27km_m300pback(108,2)/3600 16],...
            [-10000 -10000 40000 40000],[.9 .9 .9],'EdgeColor','none');
    end
    hold on
    ptot=plot(model_time_hr,Ma_pos_plot(4,:)/1000000,'k-');
    ax = gca;
    ax.ColorOrderIndex = 1;
    phfp=plot(model_time_hr,Ma_pos_plot(1,:)/1000000);
    ax = gca;
    ax.ColorOrderIndex = 5;
    pifp=plot(model_time_hr,Ma_pos_plot(2,:)/1000000);
    ax = gca;
    ax.ColorOrderIndex = 2;
    prfp=plot(model_time_hr,Ma_pos_plot(3,:)/1000000);
    ax = gca;
    ax.ColorOrderIndex = 1;
    plot(model_time_hr,Ma_neg_plot(4,:)/1000000,'k--');
    plot(model_time_hr,Ma_neg_plot(1,:)/1000000,'--');
    ax = gca;
    ax.ColorOrderIndex = 5;
    plot(model_time_hr,Ma_neg_plot(2,:)/1000000,'--');
    ax = gca;
    ax.ColorOrderIndex = 2;
    plot(model_time_hr,Ma_neg_plot(3,:)/1000000,'--');
    xlim([min_x max_x])
    %xlim([0 40])
    ylim([-3 3]);
    %set(gca,'Ytick',[0 0.025 0.05 0.075]);
    if b == 3
         ylabel({'Total';'mass flux';'(kg/s*10^6)'})
%          columnlegend(2,{'+total','+HFP','+IFP','+RFP','-total',...
%              '-HFP','-IFP','-RFP'},'Location','NorthEast');
         %xlabel('Probability');
    elseif b == 4
         %xlabel('Time (hr)');
         %set(gca,'yticklabel','');
    end
    set(gca,'Fontsize',14);
    %box(sp1,'off')
    set(gca,'xticklabel',[])
    set(gca,'Layer','top')
    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

%% Set positions of subplots 3 to 6

sp1pos = get(sp1,'Position');
sp3pos = get(sp3,'Position');
sp4pos = get(sp4,'Position');
sp5pos = get(sp5,'Position');
sp6pos = get(sp6,'Position');
sp7pos = get(sp7,'Position');
sp1pos(1) = 0.15;
sp3pos(1) = 0.15;
sp5pos(1) = 0.15;
sp7pos(1) = 0.15;
sp3pos(2) = 0.5176;
sp4pos(2) = 0.5176;
sp3pos(4) = 0.1884;
sp4pos(4) = 0.1884;
sp5pos(4) = 0.1884;
sp6pos(4) = 0.1884;
set(sp1,'Position',sp1pos);
set(sp3,'Position',sp3pos);
set(sp4,'Position',sp4pos);
set(sp5,'Position',sp5pos);
set(sp6,'Position',sp6pos);
set(sp7,'Position',sp7pos);

subplot(sp1);
set(get(get(sp1p1,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','off'); % Exclude line from legend
sp1_legend = {['Flood peak' char(10) 'in reach'],'Right FP','Left FP'};
legend(sp1_legend,'Location','east');
subplot(sp5);
set(get(get(pfi1,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','off'); % Exclude line from legend
sp5_legend = {'+ total','+ HFP','+ IFP','+ RFP'};
legend(sp5_legend,'Location','east');
subplot(sp6);
set(get(get(pfi2,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','off'); % Exclude line from legend
set(get(get(ptot,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','off'); % Exclude line from legend
set(get(get(phfp,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','off'); % Exclude line from legend
set(get(get(pifp,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','off'); % Exclude line from legend
set(get(get(prfp,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','off'); % Exclude line from legend
sp6_legend = {[char(8722) ' total'],[char(8722) ' HFP'],...
    [char(8722) ' IFP'],[char(8722) ' RFP']};
legend(sp6_legend,'Location','east');

print('-fillpage',['/Users/Colin/Documents/MATLAB/DFlow/'...
    'mass_Mo_chap3/figures/Massflux_v_time_wtype_8_5x9'],'-dpdf');