%% Reproduce Figure 4, which plots the scaled standard deviation range as a
% a function of AR(1) coefficient and the length of the reference interval.

%As the reference interval is positioned symmetrically about the midpoint
%of the time series, Delta is the difference in standard deviations between
%the midpoint and the first time point, after calculating the anomalies,
%scaled by the original standard deviation


%set sigma to one:
sig=1;

alph_vec=0:0.0125:0.975;
K_av=length(alph_vec);

N_total=252;

Midder=mean(1:1:N_total);
N_s_Vec=2:2:126;


K_NsV=length(N_s_Vec);

%% to avoid nested loops, calculate a distance matrix first, giving the separations: 
Dist_Mat=NaN(N_total);
for ii=1:1:N_total
    for jj=1:1:N_total
      Dist_Mat(ii,jj)=abs(ii-jj);
    end
end

%% run the loop. note that we only need to calculate the original cov mat
%once for each value of alpha. 

%empty matrices:
diff_matrix=NaN(K_av, K_NsV);
%difference in std from midpoint to endpoint
diff_matrix_scaled=NaN(K_av, K_NsV);
%difference in std from midpoint to endpoint, scaled by original std. 


for JJ=1:1:K_av
    
    Omega=(sig^2*alph_vec(JJ).^Dist_Mat)/(1-alph_vec(JJ)^2);

    for KK=1:1:K_NsV
        
        inds_lims=Midder+[-1, 1]*(N_s_Vec(KK)-1)/2;
        
        
        AA=zeros(N_total);
        AA(:, inds_lims(1):1:inds_lims(2))=-1/length(inds_lims(1):1:inds_lims(2));
        Omega_Anom = (eye(N_total)+AA)*Omega*(eye(N_total)+AA)';

        diff_matrix(JJ,KK)=sqrt(Omega_Anom(1,1))-sqrt(Omega_Anom(floor(Midder), floor(Midder)));
        
        diff_matrix_scaled(JJ,KK)=diff_matrix(JJ,KK)/sqrt(Omega(1,1));
        

        clear inds_lims AA Omega_Anom
        
    end
    clear Omega
    JJ
end




%% better
N_s_half=diff(N_s_Vec(1:1:2))/2;
Alph_half=diff(alph_vec(1:1:2))/2;

%set colormap
colormap('default')
ccmap=colormap;

ColLims_T=max((diff_matrix_scaled(:)))*[0, 1];

cN=length(ccmap);
colormap(ccmap)


figure(10), clf
colormap(ccmap)
fig=figure(10);
set(fig,'DoubleBuffer','on', 'color',  [1.0 1.0 1.0], 'position', [10,2,1000,790], 'renderer', 'painters');
set(gca, 'fontSize', 15)


for JJ=1:1:K_av
    for KK=1:1:K_NsV
        
        indi=max(1, ceil((diff_matrix_scaled(JJ,KK)-ColLims_T(1))/diff(ColLims_T)*cN));
        
        lon_mat=[(+N_s_half+N_s_Vec(KK))*ones(2,1);(-N_s_half+N_s_Vec(KK))*ones(2,1);]';
        
        lat_mat=[(-Alph_half+alph_vec(JJ))*ones(1,1);(+Alph_half+alph_vec(JJ))*ones(2,1);(-Alph_half+alph_vec(JJ))*ones(1,1)]';
        k1=patch(fliplr(lon_mat), fliplr(lat_mat), ccmap(indi,:), 'edgecolor', ccmap(indi,:)); hold on
        k2=plot([lon_mat, lon_mat], [lat_mat, lat_mat], 'k-', 'Color', ccmap(indi,:)); hold on

    end
end

%

contour(N_s_Vec, alph_vec, diff_matrix_scaled, 0.1923, 'color', 'k', 'linewidth', 2)

set(gca, 'Layer', 'top') 
fact=1.0;
axis([-fact*N_s_half+N_s_Vec(1), +fact*N_s_half+N_s_Vec(end), -fact*Alph_half+alph_vec(1), +fact*Alph_half+alph_vec(end)])

plot(30, 0.3430, 'w+', 'linewidth', 3, 'markersize', 8), hold on
set(gca, 'box', 'on')
    caxis('manual')
    caxis(ColLims_T)
    
pp=colorbar('peer', gca, 'vert', 'linewidth', 1, 'YAxisLocation', 'right', 'fontsize', 15);


xlabel('Length of the reference interval, N*')
ylabel('Value of the AR(1) coefficient, \alpha')

