
%% load the data
data=load('yearly_outside_temperatures.dat');
t=data(:,1);   % time in seconds
d=t/24/60/60;  % day number
d=d-min(d);    % removing large unnecessary offset
d=d+90;        % offset to match Matlab calendar  to unix one
to=data(:,2);  % outside temperature
Ndata=length(to);  % how many points

figure(1);
plot(d, to, 'x');
title('Raw data')
datetick('x','mmm'); % x to date

%% initial guessed fit
pin_simple=[7, 0, 10]; % [amplitude, phase, median_temperature]
fh1=figure(1); % get new figure handle
plot_data_and_fit(d, to, @temp_vs_day_simple, pin_simple);
title('Raw data');
datetick('x','mmm'); % x to date


%% We will not fit raw data even with simple model.
% If we attempt to fit the large data set it will take forever
% since it very time consuming calculate chi^2 
% so I commented out fit of raw data
% [pout_simple, tof, res] = fitter(d,to,@temp_vs_day_simple, pin_simple);

%% Have a look at the temperature data vs its  index
% data vector is very long several 100 thousands data points
fh2=figure(2); % get new figure handle
plot(to);
title('Raw data');

%% let's reduce the data size for now with running average
% the size of yearly variation  is about 10'000 points
% so we are save with filter size in the order of 100 points or larger
Nf=1500;
dr=reduce_with_average(d,Nf);
tor=reduce_with_average(to,Nf);
figure(1)
plot_data_and_fit(dr, tor, @temp_vs_day_simple, pin_simple);
title('Filtered data');
datetick('x','mmm'); % x to date

%% now let's fit the filtered data with simplified temperature model
%  i. e. just simple cos with offset
[pout_simple, torf, res] = fitter(dr,tor,@temp_vs_day_simple, pin_simple);

%let's transfer parameter vector to something human readable
disp('yearly temperature amplitude=');
disp(pout_simple(1));
disp('median temperature=');
disp(pout_simple(3));

% output of fitted data
figure(1)
plot_data_and_fit(dr, tor, @temp_vs_day_simple, pout_simple);
title('Filtered data');
datetick('x','mmm'); % x to date


%% may be filtered too much let's  have a look at the unfiltered data
figure(2)
[tof, resf]=plot_data_and_fit(d, to, @temp_vs_day_simple, pout_simple);
title('Raw data');
datetick('x','mmm'); % x to date

%% lets focus on residuals from RAW data but we also will take a point only every Nstep points
% otherwise our fit will take forever again

% residuals became our new data
Nstep=100;
ds=  d(1:Nstep:Ndata);
tos=to(1:Nstep:Ndata);
resfs=resf(1:Nstep:Ndata);

figure(3)
plot(ds,resfs);
title('Raw daily variations');

%% let see at the initially guessed fit
pin_daily=[8, -pi/2];
plot_data_and_fit(ds, resfs, @temp_vs_day_time, pin_daily);
title('Raw daily variations');
datetick('x','mmm'); % x to date



%% we fit data with daily variation model
% also just a simple cos but with no offset (since it is taken care by yearly variations model)
[pout_daily, torf, res] = fitter(ds,resfs,@temp_vs_day_time, pin_daily);
disp('daily temperature amplitude=');
disp(pin_daily(1));

figure(3)
plot_data_and_fit(ds, resfs, @temp_vs_day_time, pout_daily);
title('Fitted daily variations');
datetick('x','mmm'); % x to date

