%plot up the histograms of the shift results only (2D and zonal mean);
clear all; close all;


%replace with load('ASCI_results.mat') that result
load('ASCI_results_original.mat', 'DEL_P_EXT','DEL_P_ASCI', 'DEL_P_S', 'DEL_P_I', 'DEL_P_C', 'TPER_ASCI', 'ZPER_ASCI', 'TPER_ASCI2', 'ZPER_ASCI2','TPER_ASCI3', 'ZPER_ASCI3', 'do_lgm', 'do_4x', 'do_hol', 'do_hos', 'ilat', 'lon', 'SH_I', 'CON_I','INT_I', 'SH_I_T', 'CON_I_T', 'INT_I_T')

%%% PLOT HISTORAMS --- ASCI ONLY FIRST

FSz=24; FN='Arial'; 
edges=(0:2.5:100)/100;
mbin=(edges(1:end-1)+edges(2:end))/2;
ZALL=ZPER_ASCI(:,1);
ZSHIFT=ZPER_ASCI(:,2);
ZCON=ZPER_ASCI(:,3);
ZINT=ZPER_ASCI(:,4);
ZINT=min([ZINT ZALL]'); %for those models where the intensification exceeds the total
[NUM_shift, edges_shift]= histcounts(ZSHIFT, edges);
[NUM_con, edges_con]= histcounts(ZCON, edges); %as difference between total and total without contraction
[NUM_int, edges_int]= histcounts(ZINT, edges);
[NUM_all, edges_all]= histcounts(ZALL, edges);

%%%% add in some special simulations to highlight
DO=do_4x; %domain of special simulations
[SNUM_shift, edges_shift]= histcounts(ZSHIFT(do_hos), edges);
[SNUM_con, edges_con]= histcounts(ZCON(do_4x), edges); %as difference between total and total without contraction
[SNUM_int, edges_int]= histcounts(ZINT(do_lgm), edges);
[SNUM_all, edges_all]= histcounts(ZALL(DO), edges);
%[NUM_weight, edges_weight]= histcounts(WPER(:), edges);

figure(19);clf; LWE=2;
NE=length(edges);
off=0.001;
vert1=20; %above the shift
vert2=12.5; %above the contraction, and intensification
vert3 = 12.5; %above the all three
alph=0.5; Salph= 0.8;
wcol=[0.1 0.6 0.1];
SRED = [.6 .2 .2]; %color for special runs -- reddish
SPURP = [.6 .4 .6]; %special purple for contraction
SBLUE = [.2 .8 .8]; %special bluish for intensification 
for ne=1:NE-1;
hold on; plot(edges(ne)+off*[1 1], [0 NUM_shift(ne)], 'r', 'LineWidth', LWE);
hold on; plot(edges(ne)+off*[1 1], [0 SNUM_shift(ne)], 'Color',SRED, 'LineWidth', LWE);
hold on; plot(edges(ne)+off*[1 1], vert1+[0 NUM_con(ne)], 'm', 'LineWidth', LWE);
hold on; plot(edges(ne)+off*[1 1], vert1+[0 SNUM_con(ne)], 'Color',SPURP, 'LineWidth', LWE);
hold on; plot(edges(ne)+off*[1 1], (vert1+vert2)+[0 NUM_int(ne)], 'b', 'LineWidth', LWE);
hold on; plot(edges(ne)+off*[1 1], (vert1+vert2)+[0 SNUM_int(ne)], 'Color',SBLUE, 'LineWidth', LWE);
hold on; plot(edges(ne)+off*[1 1], (2*vert2+vert1)+[0 NUM_all(ne)], 'k', 'LineWidth', LWE);

%hold on; plot(edges(ne)+off*[1 1], (2*vert2+vert1+vert3)+[0 NUM_weight(ne)], 'Color', wcol, 'LineWidth', LWE);
end
for ne=2:NE;
    hold on; plot(edges(ne)-off*[1 1], (2*vert2+vert1)+[0 NUM_all(ne-1)], 'k', 'LineWidth', LWE);
    hold on; plot([edges(ne-1)+off edges(ne)-off], (2*vert2+vert1)+NUM_all(ne-1)*[1 1], 'k', 'LineWidth', LWE)
    patch([edges(ne-1)+off edges(ne-1)+off edges(ne)-off edges(ne)-off], (2*vert2+vert1)+[0 NUM_all(ne-1)*[1 1] 0], 'k', 'FaceAlpha', alph, 'EdgeColor', 'k')

    
    hold on; plot(edges(ne)-off*[1 1], (vert1+vert2)+[0 NUM_int(ne-1)], 'b', 'LineWidth', LWE);
    hold on; plot([edges(ne-1)+off edges(ne)-off], (vert1+vert2)+NUM_int(ne-1)*[1 1], 'b', 'LineWidth', LWE)
    patch([edges(ne-1)+off edges(ne-1)+off edges(ne)-off edges(ne)-off], (vert1+vert2)+[0 NUM_int(ne-1)*[1 1] 0], 'b', 'FaceAlpha', alph, 'EdgeColor', 'b')
    
    hold on; plot(edges(ne)-off*[1 1], (vert1+vert2)+[0 SNUM_int(ne-1)], 'Color', SBLUE, 'LineWidth', LWE);
    hold on; plot([edges(ne-1)+off edges(ne)-off], (vert1+vert2)+SNUM_int(ne-1)*[1 1], 'Color', SBLUE, 'LineWidth', LWE)
    patch([edges(ne-1)+off edges(ne-1)+off edges(ne)-off edges(ne)-off], (vert1+vert2)+[0 SNUM_int(ne-1)*[1 1] 0], 'b', 'FaceAlpha', Salph, 'EdgeColor', SBLUE, 'FaceColor', SBLUE)
    
    
    hold on; plot(edges(ne)-off*[1 1], vert1+[0 NUM_con(ne-1)], 'm', 'LineWidth', LWE);
    hold on; plot([edges(ne-1)+off edges(ne)-off], vert1+NUM_con(ne-1)*[1 1], 'm', 'LineWidth', LWE);
    patch([edges(ne-1)+off edges(ne-1)+off edges(ne)-off edges(ne)-off], vert1+[0 NUM_con(ne-1)*[1 1] 0], 'm', 'FaceAlpha', alph, 'EdgeColor', 'm')
    
    hold on; plot(edges(ne)-off*[1 1], vert1+[0 SNUM_con(ne-1)], 'Color',SPURP, 'LineWidth', LWE);
    hold on; plot([edges(ne-1)+off edges(ne)-off], vert1+SNUM_con(ne-1)*[1 1], 'Color',SPURP, 'LineWidth', LWE);
    patch([edges(ne-1)+off edges(ne-1)+off edges(ne)-off edges(ne)-off], vert1+[0 SNUM_con(ne-1)*[1 1] 0], 'm', 'FaceAlpha', Salph, 'EdgeColor', SPURP, 'FaceColor', SPURP)
    
    
    hold on; plot(edges(ne)-off*[1 1], [0 NUM_shift(ne-1)], 'r', 'LineWidth', LWE);
    hold on; plot([edges(ne-1)+off edges(ne)-off], NUM_shift(ne-1)*[1 1], 'r', 'LineWidth', LWE);
    patch([edges(ne-1)+off edges(ne-1)+off edges(ne)-off edges(ne)-off], [0 NUM_shift(ne-1)*[1 1] 0], 'r', 'FaceAlpha', alph, 'EdgeColor', 'r')
    
     hold on; plot(edges(ne)-off*[1 1], [0 SNUM_shift(ne-1)], 'Color',SRED, 'LineWidth', LWE);
     hold on; plot([edges(ne-1)+off edges(ne)-off], SNUM_shift(ne-1)*[1 1], 'Color',SRED, 'LineWidth', LWE);
     patch([edges(ne-1)+off edges(ne-1)+off edges(ne)-off edges(ne)-off], [0 SNUM_shift(ne-1)*[1 1] 0], 'r', 'FaceAlpha', Salph, 'EdgeColor', SRED, 'FaceColor', SRED)
%     
    
   % %seasonal weight shading etc
   % hold on; plot(edges(ne)-off*[1 1], (2*vert2+vert1+vert3)+[0 NUM_weight(ne-1)], 'Color',wcol, 'LineWidth', LWE);
   % hold on; plot([edges(ne-1)+off edges(ne)-off], (2*vert2+vert1+vert3)+NUM_weight(ne-1)*[1 1], 'Color',wcol, 'LineWidth', LWE)
   % patch([edges(ne-1)+off edges(ne-1)+off edges(ne)-off edges(ne)-off], (2*vert2+vert1+vert3)+[0 NUM_weight(ne-1)*[1 1] 0], 'g','FaceColor', wcol, 'EdgeColor', wcol, 'FaceAlpha', alph)
end
hold on; plot([min(edges) max(edges)], [0 0], 'r-', 'LineWidth', LWE) 
hold on; plot([min(edges) max(edges)], vert1+[0 0], 'm-', 'LineWidth', LWE) 
hold on; plot([min(edges) max(edges)], (vert1+vert2)+[0 0], 'b-', 'LineWidth', LWE) 
hold on; plot([min(edges) max(edges)], (vert1+2*vert2)+[0 0], 'k-', 'LineWidth', LWE)
%hold on; plot([min(edges) max(edges)], (vert1+2*vert2+vert3)+[0 0], '-', 'Color', wcol, 'LineWidth', LWE)
ymn=0;  ymx=55;%ymx=67;
xmn=min(edges); xmx=max(edges);
axis([xmn xmx ymn ymx]);
tc=[0 5 10];
ytic=[tc (vert1+tc) (vert1+vert2+tc) (vert1+2*vert2+tc)];
xtic=0:.1:1;
set(gca, 'FontSize', FSz, 'YTick', ytic, 'FontName', FN, 'LineWidth', 3, 'YTickLabel', [], 'Xtick', xtic); box on; grid on

LWM=3.5; %line width for the ensemble mean 
hold on; plot(mean(ZSHIFT)*[1 1], [ymn ymx], 'r--', 'LineWidth', LWM);
hold on; plot(mean(ZCON)*[1 1], [ymn ymx], 'm--', 'LineWidth', LWM);
hold on; plot(mean(ZINT)*[1 1], [ymn ymx], 'b--', 'LineWidth', LWM);
hold on; plot(mean(ZALL)*[1 1], [ymn ymx], 'k--', 'LineWidth', LWM);
%hold on; plot(mean(WPER)*[1 1], [ymn ymx], '--','Color', wcol, 'LineWidth', LWM);


%%%%%%%%%%%%%%%%%%%%%%%% and 2d
ALL=TPER_ASCI(:,1);
SHIFT=TPER_ASCI(:,2);
CON=TPER_ASCI(:,3);
INT=TPER_ASCI(:,4);
INT=min([INT ALL]'); %for those models where the intensification exceeds the total
[NUM_shift, edges_shift]= histcounts(SHIFT, edges);
[NUM_con, edges_con]= histcounts(CON, edges); %as difference between total and total without contraction
[NUM_int, edges_int]= histcounts(INT, edges);
[NUM_all, edges_all]= histcounts(ALL, edges);
%[NUM_weight, edges_weight]= histcounts(WPER(:), edges);

figure(20);clf; 
NE=length(edges);
off=0.001;
vert1=40; %above the shift
vert2=25; %above the contraction, and intensification
vert3 = 17.5; %above the all three
alph=0.5;
wcol=[0.1 0.6 0.1];
for ne=1:NE-1;
hold on; plot(edges(ne)+off*[1 1], [0 NUM_shift(ne)], 'r', 'LineWidth', LWE);
hold on; plot(edges(ne)+off*[1 1], vert1+[0 NUM_con(ne)], 'm', 'LineWidth', LWE);
hold on; plot(edges(ne)+off*[1 1], (vert1+vert2)+[0 NUM_int(ne)], 'b', 'LineWidth', LWE);
hold on; plot(edges(ne)+off*[1 1], (2*vert2+vert1)+[0 NUM_all(ne)], 'k', 'LineWidth', LWE);

%hold on; plot(edges(ne)+off*[1 1], (2*vert2+vert1+vert3)+[0 NUM_weight(ne)], 'Color', wcol, 'LineWidth', LWE);
end
for ne=2:NE;
    hold on; plot(edges(ne)-off*[1 1], (2*vert2+vert1)+[0 NUM_all(ne-1)], 'k', 'LineWidth', LWE);
    hold on; plot([edges(ne-1)+off edges(ne)-off], (2*vert2+vert1)+NUM_all(ne-1)*[1 1], 'k', 'LineWidth', LWE)
    patch([edges(ne-1)+off edges(ne-1)+off edges(ne)-off edges(ne)-off], (2*vert2+vert1)+[0 NUM_all(ne-1)*[1 1] 0], 'k', 'FaceAlpha', alph, 'EdgeColor', 'k')

    
    hold on; plot(edges(ne)-off*[1 1], (vert1+vert2)+[0 NUM_int(ne-1)], 'b', 'LineWidth', LWE);
    hold on; plot([edges(ne-1)+off edges(ne)-off], (vert1+vert2)+NUM_int(ne-1)*[1 1], 'b', 'LineWidth', LWE)
    patch([edges(ne-1)+off edges(ne-1)+off edges(ne)-off edges(ne)-off], (vert1+vert2)+[0 NUM_int(ne-1)*[1 1] 0], 'b', 'FaceAlpha', alph, 'EdgeColor', 'b')
    
    hold on; plot(edges(ne)-off*[1 1], vert1+[0 NUM_con(ne-1)], 'm', 'LineWidth', LWE);
    hold on; plot([edges(ne-1)+off edges(ne)-off], vert1+NUM_con(ne-1)*[1 1], 'm', 'LineWidth', LWE);
    patch([edges(ne-1)+off edges(ne-1)+off edges(ne)-off edges(ne)-off], vert1+[0 NUM_con(ne-1)*[1 1] 0], 'm', 'FaceAlpha', alph, 'EdgeColor', 'm')
    
    hold on; plot(edges(ne)-off*[1 1], [0 NUM_shift(ne-1)], 'r', 'LineWidth', LWE);
    hold on; plot([edges(ne-1)+off edges(ne)-off], NUM_shift(ne-1)*[1 1], 'r', 'LineWidth', LWE);
    patch([edges(ne-1)+off edges(ne-1)+off edges(ne)-off edges(ne)-off], [0 NUM_shift(ne-1)*[1 1] 0], 'r', 'FaceAlpha', alph, 'EdgeColor', 'r')
    
   % %seasonal weight shading etc
   % hold on; plot(edges(ne)-off*[1 1], (2*vert2+vert1+vert3)+[0 NUM_weight(ne-1)], 'Color',wcol, 'LineWidth', LWE);
   % hold on; plot([edges(ne-1)+off edges(ne)-off], (2*vert2+vert1+vert3)+NUM_weight(ne-1)*[1 1], 'Color',wcol, 'LineWidth', LWE)
   % patch([edges(ne-1)+off edges(ne-1)+off edges(ne)-off edges(ne)-off], (2*vert2+vert1+vert3)+[0 NUM_weight(ne-1)*[1 1] 0], 'g','FaceColor', wcol, 'EdgeColor', wcol, 'FaceAlpha', alph)
end
hold on; plot([min(edges) max(edges)], [0 0], 'r-', 'LineWidth', LWE) 
hold on; plot([min(edges) max(edges)], vert1+[0 0], 'm-', 'LineWidth', LWE) 
hold on; plot([min(edges) max(edges)], (vert1+vert2)+[0 0], 'b-', 'LineWidth', LWE) 
hold on; plot([min(edges) max(edges)], (vert1+2*vert2)+[0 0], 'k-', 'LineWidth', LWE)
%hold on; plot([min(edges) max(edges)], (vert1+2*vert2+vert3)+[0 0], '-', 'Color', wcol, 'LineWidth', LWE)
ymn=0;  ymx=vert1+2*vert2+vert3;%ymx=67;
xmn=min(edges); xmx=max(edges);
axis([xmn xmx ymn ymx]);
tc=[0 5 10 15 20];
ytic=[tc (vert1+tc) (vert1+vert2+tc) (vert1+2*vert2+tc)];
set(gca, 'FontSize', FSz, 'YTick', ytic, 'FontName', FN, 'LineWidth', 3, 'YTickLabel', [], 'XTick', xtic); box on; grid on

hold on; plot(mean(SHIFT)*[1 1], [ymn ymx], 'r--', 'LineWidth', LWM);
hold on; plot(mean(CON)*[1 1], [ymn ymx], 'm--', 'LineWidth', LWM);
hold on; plot(mean(INT)*[1 1], [ymn ymx], 'b--', 'LineWidth', LWM);
hold on; plot(mean(ALL)*[1 1], [ymn ymx], 'k--', 'LineWidth', LWM);
%hold on; plot(mean(WPER)*[1 1], [ymn ymx], '--','Color', wcol, 'LineWidth', LWM);