function [t,T,R,Tanom,Trms,Tady,Rndy,c,s]=seasonality9(plotit,P,nday)
% [t,T,R,Tanom,Trms,Tady,Rndy,c,s]=SEASONALITY9(plotit,P,nday)
%
% INPUT:
%
% plotit   0 calculates a correlation
%          1 also plots the data
% P        The periods in days (e.g., [181 365])
% nday     The summary window in days (e.g., 3*28)
%
% OUTPUT:
%
% t        times
% T        temperatures
% R        precipitation
% Tanom    temperature anomalies
% Trms     mean-squared prediction error
% Tady     temperature anomalies AVERAGED over n days
% Rndy     precipitation SUMMED over n days
% c        correlation matrix
% s        significance matrix
%
% EXAMPLES:
%
% seasonality9(0,[181 365],3*28);
% seasonality9(1,[181 365],3*28);
% [t,T,R,Tanom,Trms,Tady,Rndy,cc,ss]=seasonality9(0,[181 365],3*28);
% for index=1:365
%   [~,~,~,~,~,~,~,c(index),s(index)]=seasonality9(0,[181 365],index);
% end
%
% Last modified by fjsimons-at-alum.mit.edu, 11/12/2024

% Load and calculate using the previous function but do not plot it
[t,T,R,Tanom,Trms,Tady,Rndy]=seasonality8(0,P,nday);

% So now we investigate the correlation between the temperature anomaly and precipitation
% summarized as averages or sums over nonoverlapping nday windows
% Prepare for maybe offsetting this by a set number?
offset=0;
x=Tady(1:end-offset);
% Demean the data
y=Rndy(1+offset:end)-mean(Rndy);

% Calculate and report the correlation
[cc,ss]=corrcoef(x,y); c=cc(2); s=ss(2);

if plotit==1
    clf
    p3=plot(x,y,'o');
    p3.MarkerFaceColor='b';
    p3.MarkerEdgeColor='b';
    p3.MarkerSize=2;
    xlabel(sprintf(sprintf('mean temperature anomaly rms %s periods %s over %s days','%4.2f',...
                          deblank(repmat('%i ',1,length(P))),'%i'),Trms,P,nday))
    ylabel(sprintf('precipitation anomaly over %i days',nday))
    grid on
    title(sprintf('correlation coefficient %4.2f significance %4.2e',c,s))
end


