%% Calculate and plot the slope data for the three regions:


load BlockAve_Slopes_win100

ALL_slopes=BlockAve_Slopes_win100;
clear BlockAve_Slopes_win100
win_len=100;

load Inputs/BARCAST_INPUT_vHNL_Jan2012 %BARCAST_INPUT
timeline=BARCAST_INPUT.Data_timeline;

first_years=timeline(1:1:end-(win_len-1));
mid_years=first_years+(win_len-rem(win_len+1,2))/2;

last_year=2011;

to_cut=length(timeline)-find(timeline==last_year);
timeline=timeline(1:1:end-to_cut);
ALL_slopes=ALL_slopes(:, 1:1:end-to_cut);
first_years=first_years(:, 1:1:end-to_cut);
mid_years=mid_years(:, 1:1:end-to_cut);

%% get the percentiles of the slopes, through time, and plot them:
pers3=round([0.05, 0.5, 0.95]*length(ALL_slopes(:, 1)));
ALL_slopes_S=sort(ALL_slopes,1);
ALL_slopes_pers3=ALL_slopes_S(pers3, :);



%% for each ensemble member, find the midyear corresponding to the largest, smallest, and
% most extreme slope:
N_samps=length(ALL_slopes(:,1));

%largest, smallest, absolute
Ex_Slopes_ALL=NaN(N_samps, 3);

for KK=1:1:N_samps
    %all
    [val, ind]=max(ALL_slopes(KK, :));
    Ex_Slopes_ALL(KK, 1)=mid_years(ind);    
    [val, ind]=min(ALL_slopes(KK, :));
    Ex_Slopes_ALL(KK, 2)=mid_years(ind);    
    [val, ind]=max(abs(ALL_slopes(KK, :)));
    Ex_Slopes_ALL(KK, 3)=mid_years(ind);
end




%% Time series, uncertainties, and histograms. 

axer=[1445 1966 -1.65 1.65];
fig=figure(201); clf
set(fig, 'position', [10,50,1300,750]);
fill_trip=[0.7, 0.7, 0.95];
g(1)=subplot(2,1,1);
set(gca, 'fontsize', 15)

[val_L, I_L]=min(ALL_slopes_pers3(2,:));
S_Ind=I_L;
plot(mid_years(S_Ind:1:end), -ones(size(mid_years(S_Ind:1:end)))*100*ALL_slopes_pers3(1,I_L), 'r-', 'linewidth', 1.5), hold on
plot(mid_years(S_Ind:1:end), -ones(size(mid_years(S_Ind:1:end)))*100*ALL_slopes_pers3(2,I_L), 'r-', 'linewidth', 1.5), hold on
plot(mid_years(S_Ind:1:end), -ones(size(mid_years(S_Ind:1:end)))*100*ALL_slopes_pers3(3,I_L), 'r-', 'linewidth', 1.5), hold on


fill([mid_years, fliplr(mid_years)], [ALL_slopes_pers3(1,:),fliplr(ALL_slopes_pers3(3,:))]*100, fill_trip, 'linewidth', 2, 'edgecolor', fill_trip), hold on
plot(mid_years, zeros(size(mid_years)),  'k-', 'linewidth', 2.0, 'color', [0,0,1]*1.0), hold on
plot(mid_years, ALL_slopes_pers3(2,:)*100, 'k-', 'linewidth', 2.5), hold on
ylabel('^oC per Century')
%title('Land Areas, 45^oN-85^oN', 'fontsize', 20)
%text(1450, 1.25, 'Land Areas, 45^oN-85^oN', 'fontsize', 15)
plot(mid_years(I_L)*ones(1,2), -100*[ALL_slopes_pers3(1,I_L), ALL_slopes_pers3(3,I_L)], 'r-', 'linewidth', 1.5), hold on
plot(mid_years(I_L)*ones(1,2), +100*[ALL_slopes_pers3(1,I_L), ALL_slopes_pers3(3,I_L)], 'r-', 'linewidth', 1.5), hold on
plot(mid_years(I_L)*ones(1,2), +100*[ALL_slopes_pers3(1,I_L), -ALL_slopes_pers3(1,I_L)], 'r--', 'linewidth', 1.5), hold on
axis(axer)
plot(mid_years(S_Ind:1:end), -ones(size(mid_years(S_Ind:1:end)))*100*ALL_slopes_pers3(1,I_L), 'r--', 'linewidth', 1.5), hold on
plot(mid_years(S_Ind:1:end), -ones(size(mid_years(S_Ind:1:end)))*100*ALL_slopes_pers3(2,I_L), 'r--', 'linewidth', 1.5), hold on
plot(mid_years(S_Ind:1:end), -ones(size(mid_years(S_Ind:1:end)))*100*ALL_slopes_pers3(3,I_L), 'r--', 'linewidth', 1.5), hold on
plot(mid_years, ALL_slopes_pers3(2,:)*100, 'k-', 'linewidth', 2.5), hold on


%
g(2)=subplot(2,1,2);
set(gca, 'fontsize', 15)
set(gca, 'fontsize', 15)
plot(-1000*ones(1,2), -1000*ones(1,2), 'r-', 'linewidth', 2), hold on
plot(-1000*ones(1,2), -1000*ones(1,2), 'b-', 'linewidth', 2), hold on
plot(-1000*ones(1,2), -1000*ones(1,2), 'k-', 'linewidth', 8), hold on
gh=legend('Max.', 'Min.', '|Max.|');
set(gh, 'position', [0.1407 0.3165 0.0896 0.1058])

[cs, vs]=hist(Ex_Slopes_ALL(:,1), 1455:10:1955);
ins=find(cs>0);
for kk=1:1:length(ins)
        plot([vs(ins(kk))-4, vs(ins(kk))-4, vs(ins(kk))+4, vs(ins(kk))+4], [0, cs(ins(kk))/N_samps, cs(ins(kk))/N_samps, 0], 'r-', 'linewidth', 2), hold on
end

[cs, vs]=hist(Ex_Slopes_ALL(:,2), 1455:10:1955);
ins=find(cs>0);
for kk=1:1:length(ins)
        plot([vs(ins(kk))-4, vs(ins(kk))-4, vs(ins(kk))+4, vs(ins(kk))+4], [0, cs(ins(kk))/N_samps, cs(ins(kk))/N_samps, 0], 'b-', 'linewidth', 2), hold on
end

[cs, vs]=hist(Ex_Slopes_ALL(:,3), 1455:10:1955);
ins=find(cs>0);
for kk=1:1:length(ins)
        fill([vs(ins(kk))-2.0, vs(ins(kk))+2.0, vs(ins(kk))+2.0, vs(ins(kk))-2.0], [0, 0, cs(ins(kk))/N_samps, cs(ins(kk))/N_samps], 'k-', 'edgecolor', 'k'), hold on
end


plot(1448:1:1962, ones(size(1448:1:1962))*1, 'r--'), hold on
axis([1448 1968 0 1.05])
xlabel('Year')
ylabel('Fraction'), hold on

linkaxes(g, 'x');
%linkaxes(h(2:1:4), 'y');




%% AND NOW THE MAPS:
load Slopes_100win_IndsCentral_Subset_pers3
start_year_min=1547;
start_year_max=1912;

SpaceTime_3pers=Slopes_100win_IndsCentral_Subset_pers3*100; %to get to degrees per century
Locs=BARCAST_INPUT.Master_Locs(BARCAST_INPUT.Inds_Central,:);

%get the indices for plotting the slopes:

Plot_Inds=[find(timeline==start_year_min), find(timeline==start_year_max)];

% PLOT the field of slopes for the two start years, along with the uncertainty. 
Full_SpaceTime_Matrix=squeeze(SpaceTime_3pers(:, :, 2));
Full_CI_Widths_Matrix=squeeze(SpaceTime_3pers(:, :, 3))-squeeze(SpaceTime_3pers(:, :, 1));

N_Plots=length(squeeze(SpaceTime_3pers(1, :, 1)));

%
Lat_pred_range=[45, 90];
Lon_pred_range=[-180, 180];

%find the locations that will be shown that will always be NaNs:
Pos_locs=[kron([Lon_pred_range(1)+2.5:5:Lon_pred_range(2)-2.5]', ones(size([Lat_pred_range(1)+2.5:5:Lat_pred_range(2)-2.5]'))) , kron(ones(size([Lon_pred_range(1)+2.5:5:Lon_pred_range(2)-2.5]')), [Lat_pred_range(1)+2.5:5:Lat_pred_range(2)-2.5]')];
Empty_locs=setdiff(Pos_locs, Locs, 'rows');

%set colormap
load HotCoolColormap_LONG
ccmap=HotCoolColormap_LONG;

%color to plot errors larger than BigErrCUT
BigErrCol=[0.85, 0.0, 0.85];
BigErrCUT=50.5; %set it so big it doesn't matter

ColLims_T=max(max(abs(SpaceTime_3pers(:, :, 2))))*[-1, 1];
ColLims_W=[min(min(squeeze(SpaceTime_3pers(:, :, 3))-squeeze(SpaceTime_3pers(:, :, 1)))),min([BigErrCUT, max(max(squeeze(SpaceTime_3pers(:, :, 3))-squeeze(SpaceTime_3pers(:, :, 1))))])];

cN=length(ccmap);
colormap(ccmap)

figure(301), gcf
fig=figure(301); clf
set(fig, 'position', [10,50,1400,700]);

colormap(ccmap)

%fig=figure(10);
%set(fig,'DoubleBuffer','on', 'color',  [1.0 1.0 1.0], 'position', [10,10,1350,800], 'renderer', 'painters');
    
subplot(2, 2, 1)
set(gca, 'fontSize', 15)
plot(-1000, -1000, 'g^',  'color', [0.85, 0.1, 0.1], 'markersize', 5.5, 'linewidth', 1.5), hold on %'markerFaceColor',[0.3, 0.8, 0.3]), hold on
plot(-1000, -1000, 'wo', 'color', [0.70, 0.70, 0.70], 'markersize', 4, 'markerFaceColor',[0.70, 0.70, 0.70]), hold on
plot(-1000, -1000, 'g+',  'color', [0.85, 0.0, 0.85], 'markersize', 5, 'linewidth', 1.5), hold on %, 'markerFaceColor',[0.85, 0.0, 0.85]), hold on
leg=legend('Proxy', 'Inst.', 'Both', 'location', [0.4622 0.2043 0.1043 0.1471]);

for LL=1:1:N_Plots
    subplot(2, 2, LL)
    set(gca, 'fontSize', 15)
    m_proj('stereographic','long',-60,'lat',90, 'radius', 45), hold on
    m_coast('line', 'color', [0 .0 .0], 'LineWidth', 1.25), hold on
    m_grid('tickdir','in', 'xtick', [0, 90, 180, -90], 'ytick', []), hold on

    
    
    %Find the NON NaN entries for the year:
    g_inds=find(isnan(Full_SpaceTime_Matrix(:, LL))==0);    
    for kk=1:1:length(g_inds);
        indi=max(1, ceil((Full_SpaceTime_Matrix(g_inds(kk), LL)-ColLims_T(1))/diff(ColLims_T)*cN));
        lon_mat=[(+2.5+Locs(g_inds(kk),1))*ones(2,1);(-2.5+Locs(g_inds(kk),1))*ones(2,1);]';
        lat_mat=[(-2.5+Locs(g_inds(kk),2))*ones(1,1);(+2.5+Locs(g_inds(kk),2))*ones(2,1);(-2.5+Locs(g_inds(kk),2))*ones(1,1)]';
        m_patch(fliplr(lon_mat), fliplr(lat_mat), ccmap(indi,:), 'edgecolor', ccmap(indi,:)), hold on
        m_plot([lon_mat, lon_mat], [lat_mat, lat_mat], 'k-', 'Color', ccmap(indi,:)), hold on
        
        %if the 90%CI covers zero, plot a grey x:
        
        if SpaceTime_3pers(g_inds(kk), LL, 1)<0 && SpaceTime_3pers(g_inds(kk), LL, 3)>0
            m_plot(lon_mat([1,3]), lat_mat([1,3]), 'linewidth', 1, 'color', [1,1,1]*.5), hold on
            m_plot(lon_mat([2,4]), lat_mat([2,4]), 'linewidth', 1, 'color', [1,1,1]*.5), hold on
        end
        
        
    end

    %and the NaN entries
    for kk=1:1:length(Empty_locs(:,1));
        lon_mat=[(+2.5+Empty_locs((kk),1))*ones(2,1);(-2.5+Empty_locs((kk),1))*ones(2,1);]';
        lat_mat=[(-2.5+Empty_locs((kk),2))*ones(1,1);(+2.5+Empty_locs((kk),2))*ones(2,1);(-2.5+Empty_locs((kk),2))*ones(1,1)]';
        m_patch(fliplr(lon_mat), fliplr(lat_mat), [0.7, 0.7, 0.7], 'edgecolor', [0.7, 0.7, 0.7]), hold on
    end


    m_coast('line', 'color', [0 .0 .0], 'LineWidth', 1.25), hold on
    m_grid('tickdir','in', 'xtick', [0, 90, 180, -90], 'ytick', []), hold on
    hh=title([num2str(timeline(Plot_Inds(LL))), '-', num2str(timeline(Plot_Inds(LL))+99)], 'Position', [+0, +1, 15], 'fontsize', 18);
    
    if LL==1
        caxis('manual')
        caxis(ColLims_T)
        pp=colorbar('vert', 'Position', [0.8539 0.6083 0.0133 0.2946], 'fontsize', 14, 'XAxisLocation', 'bottom', 'YAxisLocation', 'right');
        tt=text(-1.2, -0.7,'Posterior Median', 'fontsize', 18, 'rotation', 90);
    end
    
    %for the uncertainty
    subplot(2, 2, 2+LL)
    set(gca, 'fontSize', 15)
    m_proj('stereographic','long',-60,'lat',90, 'radius', 45), hold on
    m_coast('line', 'color', [0 .0 .0], 'LineWidth', 1.25), hold on
    m_grid('tickdir','in', 'xtick', [0, 90, 180, -90], 'ytick', []), hold on

    %Find the NON NaN entries for the year, small uncertainty:
    g_inds=find(isnan(Full_CI_Widths_Matrix(:, LL))==0 & Full_CI_Widths_Matrix(:, LL)<BigErrCUT);    
    for kk=1:1:length(g_inds);
        indi=max(1, ceil((Full_CI_Widths_Matrix(g_inds(kk), LL)-ColLims_W(1))/diff(ColLims_W)*cN));
        lon_mat=[(+2.5+Locs(g_inds(kk),1))*ones(2,1);(-2.5+Locs(g_inds(kk),1))*ones(2,1);]';
        lat_mat=[(-2.5+Locs(g_inds(kk),2))*ones(1,1);(+2.5+Locs(g_inds(kk),2))*ones(2,1);(-2.5+Locs(g_inds(kk),2))*ones(1,1)]';
        m_patch(fliplr(lon_mat), fliplr(lat_mat), ccmap(indi,:), 'edgecolor', ccmap(indi,:)), hold on
    end

    %Find the NON NaN entries for the year, large uncertainty:
    g_inds=find(isnan(Full_CI_Widths_Matrix(:, LL))==0 & Full_CI_Widths_Matrix(:, LL)>=BigErrCUT);
    for kk=1:1:length(g_inds);
        indi=max(1, ceil((Full_CI_Widths_Matrix(g_inds(kk), LL)-ColLims_W(1))/diff(ColLims_W)*cN));
        lon_mat=[(+2.5+Locs(g_inds(kk),1))*ones(2,1);(-2.5+Locs(g_inds(kk),1))*ones(2,1);]';
        lat_mat=[(-2.5+Locs(g_inds(kk),2))*ones(1,1);(+2.5+Locs(g_inds(kk),2))*ones(2,1);(-2.5+Locs(g_inds(kk),2))*ones(1,1)]';
        m_patch(fliplr(lon_mat), fliplr(lat_mat), [0.4, 0.7, 0.4], 'edgecolor', [0.4, 0.7, 0.4]), hold on
    end


    %and the always NaN locations
    for kk=1:1:length(Empty_locs(:,1));
        lon_mat=[(+2.5+Empty_locs((kk),1))*ones(2,1);(-2.5+Empty_locs((kk),1))*ones(2,1);]';
        lat_mat=[(-2.5+Empty_locs((kk),2))*ones(1,1);(+2.5+Empty_locs((kk),2))*ones(2,1);(-2.5+Empty_locs((kk),2))*ones(1,1)]';
        m_patch(fliplr(lon_mat), fliplr(lat_mat), [0.7, 0.7, 0.7], 'edgecolor', [0.7, 0.7, 0.7]), hold on
    end
    
        m_coast('line', 'color', [0 .0 .0], 'LineWidth', 1.25), hold on
    m_grid('tickdir','in', 'xtick', [0, 90, 180, -90], 'ytick', []), hold on
    
    %find where we have data 
    I_g_I=find(isnan(BARCAST_INPUT.Inst_Data(:, Plot_Inds(LL)))==0);
    I_g_P1=find(isnan(BARCAST_INPUT.Prox_Data1(:, Plot_Inds(LL)))==0);
    I_g_P2=find(isnan(BARCAST_INPUT.Prox_Data2(:, Plot_Inds(LL)))==0);
    I_g_P3=find(isnan(BARCAST_INPUT.Prox_Data3(:, Plot_Inds(LL)))==0);
    I_g_P=unique([I_g_P1;I_g_P2;I_g_P3]); 
    
     %only the Inst
     [Only_I, inds_Only_I]=setdiff(I_g_I, I_g_P);
     %only the Prox
     [Only_P, inds_Only_P]=setdiff(I_g_P, I_g_I);
     %Both
     [Both, inds_Both]=intersect(I_g_P, I_g_I);
    
     %plot
     %only the Prox
     m_plot(BARCAST_INPUT.Master_Locs(Only_P,1), BARCAST_INPUT.Master_Locs(Only_P,2), 'g^',  'color', [0.85, 0.1, 0.1], 'markersize', 5.5, 'linewidth', 1.5), hold on %'markerFaceColor',[0.3, 0.8, 0.3]), hold on
    %only the Inst
     m_plot(BARCAST_INPUT.Master_Locs(Only_I,1), BARCAST_INPUT.Master_Locs(Only_I,2), 'wo', 'color', [0.70, 0.70, 0.70], 'markersize', 4, 'markerFaceColor',[0.70, 0.70, 0.70]), hold on
    %Both
     m_plot(BARCAST_INPUT.Master_Locs(Both,1), BARCAST_INPUT.Master_Locs(Both,2), 'g+',  'color', [0.85, 0.0, 0.85], 'markersize', 5, 'linewidth', 1.5), hold on %, 'markerFaceColor',[0.85, 0.0, 0.85]), hold on
        
    hh=title([num2str(timeline(Plot_Inds(LL))), '-', num2str(timeline(Plot_Inds(LL))+99)], 'Position', [+0, +1, 15], 'fontsize', 18);


     if LL==1
        caxis('manual')
        caxis(ColLims_W)
        pp=colorbar('vert', 'Position', [0.8539 0.1157 0.0133 0.3615], 'fontsize', 14, 'XAxisLocation', 'bottom', 'YAxisLocation', 'right');
        tt=text(-1.2, -0.85,'90% Credible Intervals', 'fontsize', 18, 'rotation', 90);
     end
    
end
%0.8539 0.6083 0.0133 0.2946
%0.8539 0.1157 0.0362 0.3615


