function [t,T,R,Tanom,Trms]=seasonality7(plotit,P)
% [t,T,R,Tanom,Trms]=SEASONALITY7(plotit,P)
%
% INPUT:
%
% plotit   0 calculates a temperature anomaly
%          1 also plots the data
% P        The periods in days (e.g., [181 365])
%
% OUTPUT:
%
% t        times
% T        temperatures
% R        precipitation
% Tanom    temperature anomalies
% Trms     mean-squared prediction error
%
% EXAMPLES:
%
% seasonality7(0,[181 365]);
% seasonality7(1,[181 365]);
% [t,T,R,Tanom,Trms]=seasonality7(0,[181 365]);
%
% Last modified by fjsimons-at-alum.mit.edu, 11/12/2024

% Load the data from the previous function but do not plot it
[t,T,R]=seasonality6(0);

% Just number the days sequentially
td=1:length(t);
% Prepare for axis labels in decades
xtix=[1:10*365-1:length(td)];
xtil=[0:10:floor(length(td)/365)];

% Let us make the temperature anomaly by taking out the likely cycles
% using the P that was input
% Use a function! Feed it the running days, the temperature, the periods
[Tanom,Trms]=anomaly(td,T,P);

if plotit==1
    % Make a quick plot of temperature anomaly and precipitation
    clf
    s1=subplot(211);
    p1=plot(td,Tanom); axis tight; ylabel('temperature anomaly')
    xlabs(s1,xtix,xtil,0); xlabel(sprintf('years since %i',year(t(1))))
    title(sprintf(sprintf('rms %s periods %s','%4.2f',...
                          deblank(repmat('%i ',1,length(P)))),Trms,P))
    s2=subplot(212); 
    p2=plot(td,R); axis tight; ylabel('precipitation')
    xlabs(s1,xtix,xtil,0); xlabel(sprintf('years since %i',year(t(1))))
end

% Here we make a (sub)function, using all "internal" variables
function [vanom,vrms]=anomaly(t,v,p)
argm=[1./p(:)*t(:)']';
F=[ones(length(t),1) sin(2*pi*argm) cos(2*pi*argm)];
A=pinv(F)*v;
vpred=F*A;
% Now here is the temperature anomaly
vanom=v-vpred;
% Check the "misfit" or the mean-squared prediction error
vrms=sqrt(sum(vanom.^2)/length(v));

% Here is another (sub)function, simply to clean up the labels
function xlabs(ax,xtix,xtil,ex)
axis(ax); xticks(xtix); xticklabels(xtil); s1.XAxis.Exponent=ex;
