function [t,T,R,Tanom,Trms,Tady,Rndy]=seasonality8(plotit,P,nday)
% [t,T,R,Tanom,Trms,Tady,Rndy]=SEASONALITY8(plotit,P,nday)
%
% INPUT:
%
% plotit   0 calculates a temperature anomaly and the summary windows
%          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
%
% EXAMPLES:
%
% seasonality8(0,[181 365],3*28);
% seasonality8(1,[181 365],3*28);
% [t,T,R,Tanom,Trms,Tady,Rndy]=seasonality8(0,[181 365],3*28);
%
% 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]=seasonality7(0,P);

% Let us summarize the "data" in n-day non-overlapping increments
% using the nday that was input
% How many times can we do this?
nfold=floor(length(R)/nday);
remdr=rem(length(R),nday);
% Simply rearrange and SUM across the nday columns
Rndy=sum(reshape(R(1:end-remdr),nday,[]),1);
% And for the temperature, how about we look at the MEAN ANOMALY across the n days
Tady=sum(reshape(Tanom(1:end-remdr),nday,[]),1)/nday;
% And now this will be the time axis in increments of n days for plottin
tnd=1:length(Rndy);

if plotit==1
    % Make a new plot of temperature anomaly and precipitation over those n days
    clf
    s1=subplot(211);
    p1=plot(tnd,Tady); axis tight; ylabel('temperature anomaly'); grid on
    title(sprintf(sprintf('rms %s periods %s sum/average window %s days','%4.2f',...
                          deblank(repmat('%i ',1,length(P))),'%i'),Trms,P,nday))
    s2=subplot(212); 
    p2=plot(tnd,Rndy); axis tight; ylabel('precipitation'); grid on
end
