% 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);
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))

% INVERSION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now do the formal analysis
P=24;
argm=tsel/P;
% Write the model on the blackboard as T=F*A;
F=[ones(length(tsel),1) cos(2*pi*argm) sin(2*pi*argm)];
 
% And this here is how you "invert" the relationship to come up with the
% weighting constants that are the amplitudes for each of the terms, and by
% 'applying' the inverse to the data vector - a 'linear regression'!
abc=pinv(F)*Tsel;
% Use the old nomenclature
a=abc(1);
b=abc(2);
c=abc(3);
% So now this is the prediction
Tpred=a+b*sin(2*pi*tsel/P)+c*cos(2*pi*tsel/P);
% And plot that instead of what you had before
axes(s2)
p2=plot(tsel,Tpred,'r'); 
ylabel('modeled T');
xlabel('hours')
grid on
axis(zax)
% Overlay
axes(s1); hold on
p3=plot(tsel,Tpred,'r'); 
hold off
% Calculate the residuals
resd=Tsel-Tpred;
% Calculate the root-mean-squared error
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))

