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

% MANUAL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First, make a plot
clf
s1=subplot(211);
% Plot time against temperature in hours since the first hour
p1=plot(tsel,Tsel);
ylabel('observed T');
xlabel('hours')
grid on
axis tight
% Grab whatever 'tight' means
zax=axis;

% The model is that T=a+b*sin(2*pi*t/P)+c*cos(2*pi*t/P);
s2=subplot(212);
% Make a wild guess
% Get ready to inquire about a 24-hour period in these HOURLY data
P=24;
a=28;
b=4.5;
c=0;
% So this is your prediction
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
% Apply whatever 'tight' meant
axis(zax)
% Calculate the residuals
resd=Tsel-Tpred;
% Calculate the root-mean-squared error
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))

% Play with these numbers - can you make a function that takes
% tsel, P1, [a b c] and play with the numbers?
% if you get these coefficients right, the time series will look the same

