clear all; clc;
load 'STOCKHOLMSDATA.mat' Tm;


A = zeros(length(Tm),3);
A(:,1) = 1;
A(:,2) = 1:length(Tm);
Ts = 240;
modellvard = zeros(length(Tm),1);
Omega = 2*pi/(365*24);
x = zeros(3,1);
%Ber??kna konstanterna
tic
for i = 1:length(Tm)
    A(i,3) = sin(Omega*(i - Ts));
end
toc
x = A\Tm;
%Ber??kna normen
modellvard(1:length(Tm)) = anpass(x,1:length(Tm),Omega,Ts);
resnorm = norm(Tm - modellvard);

Ts2 = 182*24;

resnorm2 = zeros(Ts2,1);
%Ber??kna normen av residualen f??r olika Ts.

for c = 1:Ts2
    %Ber??kna konstanterna med varierande Ts.
    for i = 1:length(Tm)
        A(i,3) = sin(Omega*(i - c));
    end
    x2 = A\Tm;
    modellvard(1:length(Tm)) = anpass(x2,1:length(Tm),Omega,c);
    resnorm2(c) = norm(Tm - modellvard);
end
%Hitta Ts med minst norm
[minnorm,index]=min(resnorm2);

for i = 1:length(Tm)
        A(i,3) = sin(Omega*(i - index));
end
x3 = A\Tm;

%Plotta figurer
figure(1);
%v??rden fr??n start
plot(1:length(Tm),anpass(x,1:length(Tm),Omega,Ts));
hold on; grid on;
xlabel('Timmar');
ylabel('Temperatur');
%Data
plot(1:length(Tm),Tm);
%Anpassade v??rden med minst norm
plot(1:length(Tm),anpass(x3,1:length(Tm),Omega,index),'k');
%Fasf??rskjutning i dagar
fasfor = round(index/24);
%Kallaste dagen p?? ??ret blir 26 Januari, Varmaste 27 Juli. Det sker en
%temperatur??kning med ungef??r 1 grad.

function s = anpass(x,t,Omega,Ts)
    s = x(1) + x(2).*t + x(3).*sin(Omega.*(t-Ts));
end