clear all; close all; clc;
CDHE_Data = xlsread('CDHE_Events_All.xlsx');
St_Name = ["Great Plains South","Midwest-Great Plains North",...
    "Northeast","Southeast","Southwest","Northwest"];
St_Number = [0,6,17,23,37,47,50];

for kk = 1 : length(St_Number)-1
    close all
    clear index_zone
    index_zone1 =  find(ismember(CDHE_Data(:,3),...
        (St_Number(kk)+1:St_Number(kk+1))));
    figure(kk)
    box1 = boxplot([CDHE_Data(index_zone1,1);110],...
        [CDHE_Data(index_zone1,2);4],'Colors','bcmr','OutlierSize',4,...
        'Orientation','horizontal');
    title('')
    
    hhh = findobj(gcf,'tag','Outliers');
    hhh(1,1).MarkerEdgeColor = 'r';
    hhh(2,1).MarkerEdgeColor = 'm';
    hhh(3,1).MarkerEdgeColor = 'c';
    hhh(4,1).MarkerEdgeColor = 'b';
    set(findobj(gca,'type','line'),'linew',0.9)
    set(gca,'fontsize',12)
    
    if ~exist([pwd '\Results'], 'dir')
        mkdir([pwd '\Results'])
    end
    
    xlim([0 62])
    set(gca,'xtick',[0,20,40,60],'xticklabel',[0,20,40,60])
    box off
    
    x0 = 10;
    y0 = 1;
    width = 2;
    height= 2;
    set(gcf,'units','inches','position',[x0,y0,width,height])
    
    if ~exist([pwd '\Results\Boxplots'], 'dir')
        mkdir([pwd '\Results\Boxplots'])
    end
    fpath = [pwd '\Results\Boxplots'];
    print(figure(kk),fullfile(fpath , [char(St_Name(kk)) '.jpeg']), '-djpeg','-r600')
    
    %ECDF plots
    clear indexzone2
    indexzone2 = find(ismember(CDHE_Data(:,3),...
        (St_Number(kk)+1:St_Number(kk+1))) & CDHE_Data(:,2) == 4);
    figure(kk+6)
    [ff,xx] = ecdf(CDHE_Data(indexzone2,1));
    hold on
    Future_CDf =  plot(xx,ff,'r','linewidth',1.5);
    clear indexzone2 ff xx
    indexzone2 = find(ismember(CDHE_Data(:,3),...
        (St_Number(kk)+1:St_Number(kk+1))) & CDHE_Data(:,2) == 1);
    
    [ff,xx] = ecdf(CDHE_Data(indexzone2,1));
    hold on
    Current_CDf = plot(xx,ff,'b','linewidth',1.5);
    box on
    xlim([0 60])
    
    grid off
    title('')
    xlabel('Duration[Day]')
    ylabel('')
    
    x0 = 1;
    y0 = 1;
    width = 3;
    height= 2;
    set(gcf,'units','inches','position',[x0,y0,width,height])
    check1 = 13;
    
    if ~exist([pwd '\Results\ECDFs'], 'dir')
        mkdir([pwd '\Results\ECDFs'])
    end
    fpath = [pwd '\Results\ECDFs'];
    print(figure(kk+6),fullfile(fpath , [char(St_Name(kk)) '.jpeg']), '-djpeg','-r600')
    
    % Finding regional percentiles
    for SC = 1 : 4
        clear index_zone CDHE_events
        index_zone1 =  find(ismember(CDHE_Data(:,3),(St_Number(kk)+1:St_Number(kk+1))) & CDHE_Data(:,2) == SC);
        No_Station = length(St_Number(kk)+1:St_Number(kk+1));
        CDHE_events = CDHE_Data(index_zone1,1);
        Prct_99(kk,SC) = prctile(CDHE_events,99);
        Prct_95(kk,SC) = prctile(CDHE_events,95);
        Prct_50(kk,SC) = prctile(CDHE_events,50);
    end
    
    % Finding the total number of CDHE days in US under diffrent scenarios
    for SC = 1 : 4
        clear index_Scenario index_zone1
        %        index_Scenario =  find(CDHE_Data(:,2) == SC);
        index_Scenario =  find(ismember(CDHE_Data(:,3),...
            (St_Number(kk)+1:St_Number(kk+1))) & CDHE_Data(:,2) == SC);
        CDHE_Days(kk,SC) =  sum(CDHE_Data(index_Scenario,1))./(No_Station*10);
        try
            CDHE_Days(kk,3) =  sum(CDHE_Data(index_Scenario,1))./(No_Station*check1);
        catch
        end
    end
end

% Finding the total number of CDHE day in each city under diffrent scenario
for kk = 1 : length(unique(CDHE_Data(:,3)))
    for SC = 1 : 4
        clear index_Scenario
        index_Scenario =  find(CDHE_Data(:,3) == kk & CDHE_Data(:,2) == SC);
        CDHE_Days_city(kk,SC) =  sum(CDHE_Data(index_Scenario,1))./10;
    end
end

%% Barplots (Regional)
for kk = 1 : 6
    clear data1
    data1 = CDHE_Days_city((St_Number(kk)+1:St_Number(kk+1)),:);
    upper1(kk,:) = std(data1);
    lower1(kk,:) = -std(data1);
end

for kk = 1 : 6
    figure(kk)
    figure(kk)
    bar_daily = bar(1:4,CDHE_Days(kk,:),'FaceColor','flat');
    ylim([0 79])
    bar_daily.CData(1,:) = [0 0 1];
    bar_daily.CData(2,:) = [0 1 1];
    bar_daily.CData(3,:) = [1 0 1];
    bar_daily.CData(4,:) = [1 0 0];
    hold on
    er = errorbar(1:4,CDHE_Days(kk,:),lower1(kk,1:end),upper1(kk,1:end));
    er.Color = [0 0 0];
    er.LineStyle = 'none';
    x0 = 1;
    y0 = 1;
    width = 2;
    height= 2.2;
    set(gcf,'units','inches','position',[x0,y0,width,height])
    if ~exist([pwd '\Results\Barplots'], 'dir')
        mkdir([pwd '\Results\Barplots'])
    end
    fpath = [pwd '\Results\Barplots'];
    print(figure(kk),fullfile(fpath , [char(St_Name(kk)) , '_DailyCDHE.jpeg']), '-djpeg','-r600')
end