%% plot a temperature series at a specific location.  

load Temperature_Array_3pers %Temperature_Array_3pers
load SM_Temperature_Array_3pers %_SMTemperature_Array_3pers
load Parameters_pers3 %Parameters_pers3

load Inputs/BARCAST_INPUT_vHNL_Jan2012 %BARCAST_INPUT
timeline=[BARCAST_INPUT.Data_timeline(1)-1, BARCAST_INPUT.Data_timeline];
Data_timeline=BARCAST_INPUT.Data_timeline;

Temperature_Array_3pers=Temperature_Array_3pers(:, 2:1:end, :);
SM_Temperature_Array_3pers=SM_Temperature_Array_3pers(:, 2:1:end, :);


%% get the mean/variance for each time series (without spatial
% considerations). 
%1 = alpha, the AR(1) coefficient
%2 = mu, the constant par in the linear mean of the AR(1) process
%3 = sigma2, the partial sill in the spatial covariance matrix
%4 = phi, the range parameter in the spatial covariance matrix

alph=Parameters_pers3(1,2);
mm=Parameters_pers3(2,2);
sig2=Parameters_pers3(3,2);
vv=sig2/(1-alph^2);

pers_of_proc=[.05, .5, .95];
proc_prior_3pers=norminv(pers_of_proc, mm, sqrt(vv));

% do the same for the nine-point smoothing. 
Nsm=9;
ww=hanning(Nsm)/sum(hanning(Nsm));

% get the cov matrix:
cc=NaN(Nsm,Nsm);
for ii=1:1:Nsm
    for jj=1:1:Nsm
        cc(ii,jj) = vv * alph^(abs(ii-jj));
    end
end

mm_sm=mm*ones(1,Nsm)*ww; %check: mean should not be affected
vv_sm=ww'*cc*ww;
proc_prior_SM_3pers=norminv(pers_of_proc, mm_sm, sqrt(vv_sm));


%% Plot two time series

% %find the New Brunswick time series
inds1=find(BARCAST_INPUT.Master_Locs(:,1)==-67.5 & BARCAST_INPUT.Master_Locs(:,2)==47.5); 
% middle of greenland
inds2=find(BARCAST_INPUT.Master_Locs(:,1)==-42.5 & BARCAST_INPUT.Master_Locs(:,2)==77.5); 
%North of Moscow
inds3=find(BARCAST_INPUT.Master_Locs(:,1)==37.5 & BARCAST_INPUT.Master_Locs(:,2)==57.5); 
%Southern part of europe
inds4=find(BARCAST_INPUT.Master_Locs(:,1)==7.5 & BARCAST_INPUT.Master_Locs(:,2)==47.5); 
%Stockholm
inds5=find(BARCAST_INPUT.Master_Locs(:,1)==17.5 & BARCAST_INPUT.Master_Locs(:,2)==57.5); 
%or choose your own
 
 
inds=[inds1, inds2];

N_pts=length(inds);
fill_trip=[0.7, 0.7, 0.95];
%fill_trip=[1,1,1]*.7;

fig=figure(12); clf
set(fig, 'position', [10,50,1300,750], 'renderer', 'painters');
for kk=1:1:N_pts
    %raw
    subplot(2*N_pts, 1, 2*kk-1)
    set(gca, 'fontsize', 15)
    fill([Data_timeline, fliplr(Data_timeline)], [squeeze(Temperature_Array_3pers(inds(kk), :,1)), fliplr(squeeze(Temperature_Array_3pers(inds(kk), :,3)))], fill_trip, 'linewidth', 2, 'edgecolor', fill_trip), hold on 
    
    plot(Data_timeline, zeros(size(Data_timeline)),  'k-', 'linewidth', 1.0, 'color', [0,0,1]*1.0), hold on
    plot(Data_timeline, ones(size(Data_timeline))*proc_prior_3pers(1),  'r--', 'linewidth', 1.25), hold on
    plot(Data_timeline, ones(size(Data_timeline))*proc_prior_3pers(3),  'r--', 'linewidth', 1.25), hold on

    plot(Data_timeline, squeeze(Temperature_Array_3pers(inds(kk), :,2)),  'k-', 'linewidth', 2.5), hold on
    axis([1395 2015 -4.0 4.0])
    ylabel('^oC')
    set(gca, 'xticklabel', {})

    %smoothed
    subplot(2*N_pts, 1, 2*kk)
    set(gca, 'fontsize', 15)
    fill([Data_timeline, fliplr(Data_timeline)], [squeeze(SM_Temperature_Array_3pers(inds(kk), :,1)), fliplr(squeeze(SM_Temperature_Array_3pers(inds(kk), :,3)))], fill_trip, 'linewidth', 2, 'edgecolor', fill_trip), hold on 
    plot(Data_timeline, squeeze(SM_Temperature_Array_3pers(inds(kk), :,2)),  'k-', 'linewidth', 2.5), hold on
    plot(Data_timeline, ones(size(Data_timeline))*proc_prior_SM_3pers(1),  'r--', 'linewidth', 1.25), hold on

    if kk==N_pts
        leg=legend('90%CI', 'Posterior Median', 'Process-level 90%CI','orientation', 'hor', 'location', [0.3941 0.4803 0.2354 0.0500]);
    end
    plot(Data_timeline, zeros(size(Data_timeline)),  'k-', 'linewidth', 1.0, 'color', [0,0,1]*1), hold on
    plot(Data_timeline, ones(size(Data_timeline))*proc_prior_SM_3pers(1),  'r--', 'linewidth', 1.25), hold on
    plot(Data_timeline, ones(size(Data_timeline))*proc_prior_SM_3pers(3),  'r--', 'linewidth', 1.25), hold on

    plot(Data_timeline, squeeze(SM_Temperature_Array_3pers(inds(kk), :,2)),  'k-', 'linewidth', 2.5), hold on
    axis([1395 2015 -4.0 4.0])
    set(gca, 'YAxisLocation', 'left')
    ylabel('^oC')
    text(1405, -2.8, 'Smoothed', 'fontsize', 18)
        title(['[', num2str(BARCAST_INPUT.Master_Locs(inds(kk),1)), 'E, ', num2str(BARCAST_INPUT.Master_Locs(inds(kk),2)), 'N]'], 'fontsize', 18)
    if kk==N_pts
        xlabel('Year')
    else
        set(gca, 'xticklabel', {})
    end
end
subplot(4,1,1)
text(1400, 2.8, '(a)', 'fontsize', 15)
subplot(4,1,2)
text(1400, 2.8, '(b)', 'fontsize', 15)
subplot(4,1,3)
text(1400, 2.8, '(c)', 'fontsize', 15)
subplot(4,1,4)
text(1400, 2.8, '(d)', 'fontsize', 15)

