% SAME %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rome=[41.905 12.474];
fnam='data_stream-oper.nc';
[T,t,closestnode]=readVariableFromNetCDF(fnam,'t2m',rome);
sel=sum(day(t)==[1:7],2) & month(t)==8 & year(t)==2024;
tsel=t(sel);
tsel=hours(tsel-tsel(1));
Tsel=T(sel)-273;
clf
s1=subplot(211);
p1=plot(tsel,Tsel);
ylabel('observed T');
xlabel('hours')
grid on
axis tight
zax=axis;
s2=subplot(212);
P=24;
a=28;
b=4.5;
c=0;
Tpred=a+b*sin(2*pi*tsel/P)+c*cos(2*pi*tsel/P);
% And plot that
p2=plot(tsel,Tpred,'r'); 
ylabel('modeled T');
xlabel('hours')
grid on
axis(zax)
resd=Tsel-Tpred;
rmser=sqrt(sum(resd.^2)/length(resd));
title(sprintf('a= %4.2f b= %4.2f c= %4.2f rmse= %4.2f',a,b,c,rmser))
P=24;
argm=tsel/P;
F=[ones(length(tsel),1) cos(2*pi*argm) sin(2*pi*argm)];
abc=pinv(F)*Tsel;
a=abc(1);
b=abc(2);
c=abc(3);
Tpred=a+b*sin(2*pi*tsel/P)+c*cos(2*pi*tsel/P);
axes(s2)
p2=plot(tsel,Tpred,'r'); 
ylabel('modeled T');
xlabel('hours')
grid on
axis(zax)
axes(s1); hold on
p3=plot(tsel,Tpred,'r'); 
hold off
resd=Tsel-Tpred;
rmser=sqrt(sum(resd.^2)/length(resd));
axes(s2)
title(sprintf('a= %4.2f b= %4.2f c= %4.2f rmse= %4.2f',abc(1),abc(2),abc(3),rmser))
P=[6 12 24 48];
argm=[1./P(:)*tsel(:)']';
F=[ones(length(tsel),1) cos(2*pi*argm) sin(2*pi*argm)];
A=pinv(F)*Tsel;
Tpred=F*A;
axes(s2)
p2=plot(tsel,Tpred,'r'); 
ylabel('modeled T');
xlabel('hours')
grid on
axis(zax)
axes(s1); hold on
delete(p3)
p3=plot(tsel,Tpred,'r'); 
hold off
resd=Tsel-Tpred;
rmser=sqrt(sum(resd.^2)/length(resd));
axes(s2)
title(sprintf('rmse= %4.2f',rmser))
axes(s2)
plot(tsel,Tsel-Tpred)
ylabel('T anomaly, 2 sigma');
xlabel('hours')
grid on
% Make a meaningful axis
axis([zax([1 2]) 0+[-1 1]*[zax(4)-zax(3)]])
% You could add ticks at the standard deviations of the anomaly!
s2.YTick=round(10*[mean(resd)+[-2 0 2]*std(resd)])/10;
ylim([-6 6]*std(resd))

% SEASONAL ANOMALY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Maybe go with daily? Weekly?

% Select all the data
% Select the times and refer them to the first hour
tsel=t;
tsel=hours(tsel-tsel(1));
% Select the temperatures and convert them to degrees Celsius
Tsel=T-273;

clf
s1=subplot(211);
p1=plot(tsel,Tsel);
ylabel('observed T');
xlabel('hours')
grid on
axis tight
s2=subplot(212);

% Now repeat all of this and find seasonal cycles on top of this
P=[6 12 24 48 24*181.0 24*365.0];

% The rest is almost indentical, so no more comments
argm=[1./P(:)*tsel(:)']';
F=[ones(length(tsel),1) cos(2*pi*argm) sin(2*pi*argm)];
A=pinv(F)*Tsel;
Tpred=F*A;
axes(s1); hold on
p3=plot(tsel,Tpred,'r'); 
hold off
axis tight
zax=axis;
resd=Tsel-Tpred;
rmser=sqrt(sum(resd.^2)/length(resd));
axes(s2)
plot(tsel,Tsel-Tpred)
ylabel('T anomaly, 2 sigma');
xlabel('hours')
grid on
axis([zax([1 2]) 0+[-1 1]*[zax(4)-zax(3)]])
s2.YTick=round(10*[mean(resd)+[-2 0 2]*std(resd)])/10;
ylim([-6 6]*std(resd))

