
%% 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-d(1)+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

%% yearly temperature model
temp_vs_day_simple_model=fittype( @(a, ph_y, Taver,x) -a*cos(2*pi/365*x + ph_y) + Taver);
temp_vs_day_simple=@(x,p) -p(1)*cos(2*pi/365*x + p(2)) + p(3);

%% 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 we need to reduce number of points

%% 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
% raw data set has in order of a million points we need to make it reasonably sized
Nf=floor(24*60*60/data_points_spacing); % designed to be about one averaged 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
tic; % reset the calculation timer
[fo_temp_vs_day_simple, gof] = fit(dr,tor,temp_vs_day_simple_model, 'StartPoint', pin_simple);
toc % output the elapsed time

%let's transfer parameter vector to something human readable
disp('yearly temperature amplitude=');
disp(fo_temp_vs_day_simple.a);
disp('median temperature=');
disp(fo_temp_vs_day_simple.Taver);

fitted_yearly_params=[ fo_temp_vs_day_simple.a, fo_temp_vs_day_simple.ph_y, fo_temp_vs_day_simple.Taver];

% output of fitted data
figure(1)
[tof, resf, plot_handles]=plot_weather_data_and_fit(dr, tor, temp_vs_day_simple, fitted_yearly_params);
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, fitted_yearly_params);
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 is  our new data and focus of interest
Nstep=floor(2*60*60/data_points_spacing); % we will use take a point every 2 hours
% note that I am not averaging since it will smooth temperature variations
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

% we will model it as cos as well
temp_vs_day_time_model=fittype( @(a,ph,x) -a*cos(2*pi/1*x - ph) );
temp_vs_day_time = @(x,p) temp_vs_day_time_model(p(1), p(2), x);

%% let see at the initially guessed fit
pin_daily=[8, pi/2]; % daily temperature variation amplitude and phase shift
[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 temperature offset 
% since it is taken care by yearly variations model
tic;
[fo_temp_vs_day_time, gof] = fit(ds,resfs,temp_vs_day_time_model, 'StartPoint',pin_daily);
toc
disp('daily temperature amplitude=');
disp(fo_temp_vs_day_time.a);

fitted_day_params=[fo_temp_vs_day_time.a, fo_temp_vs_day_time.ph];

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

% clearly the fit is quite bad, since sinus is not the best looking model

%% 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(fo_temp_vs_day_time.ph*24/(2*pi));

