
%% load the data
data=load('yearly_outside_temperatures.dat');
t=data(:,1);   % time in seconds (unix date format)
d=t/24/60/60;  % time conversion to day number convention
date_offset=datenum([1970 1 1 0 0 t(1)]); % unix epoch to Matlab
d=d-min(d);    % removing large unnecessary offset
d=d+date_offset; % converting to Matlab date format
to=data(:,2);  % outside temperature
Ndata=length(to);  % how many points

figure(1);
plot(d, to, '.');
ylabel('Temp (^oC)');
title('Raw data');
datetick('x','yyyy/mmm'); % x to date

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


%% Have a look at the temperature data array length 
data_length=length(to)

% 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 even to plot so many points.
% Calculation of  chi^2 will kill us
% so I commented out fit of raw data
% [pout_simple, tof, res] = fitter(d,to,@temp_vs_day_simple, pin_simple);

%% Average spacing in seconds between consequent points excluding longer then hour voids
dt=diff(t); 
data_points_spacing=mean(dt(dt <= 3600))

%% 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=60*24; % designed to be one data point per day
dr=reduce_with_average(d,Nf);
tor=reduce_with_average(to,Nf);
figure(1)
plot_weather_data_and_fit(dr, tor, @temp_vs_day_simple, pin_simple);
title('Filtered data');

%% 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)
[tof, resf, plot_handles]=plot_weather_data_and_fit(dr, tor, @temp_vs_day_simple, pout_simple);
title('Filtered data');


%% may be we filtered too much let's  have a look at the unfiltered data
figure(2)
[tof, res, plot_handles]=plot_weather_data_and_fit(d, to, @temp_vs_day_simple, pout_simple);
title('Raw data');

%% 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; % we still use every 100th point instead of all the data
ds=  d(1:Nstep:Ndata);
tos=to(1:Nstep:Ndata);
resfs=res(1:Nstep:Ndata);

figure(3)
plot(ds,resfs);
title('Raw daily variations');
datetick('x','yyyy/mmm'); % x to date

%% let see at the initially guessed fit
pin_daily=[8, pi/2];
[tof, resf, plot_handles]=plot_weather_data_and_fit(ds, resfs, @temp_vs_day_time, pin_daily);
ylabel('Temp (^oC)');
title('Raw daily variations');



%% 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(pout_daily(1));

figure(3)
[tof, resf, plot_handles]=plot_weather_data_and_fit(ds, resfs, @temp_vs_day_time, pout_daily);
title('Fitted daily variations');

% clearly the fit is quite bad

%% But we can extract a useful information from it
% looks like our fit is in sync with daily variation even though it often misses.
% So we can find the coldest hour of the day
disp('Coldest hour of the day=');
disp(pout_daily(2)*24/(2*pi));

