%uloha 1
load list_of_stars.dat
x = list_of_stars(:,1);     %hodnoty poradia = x-ova os
y = list_of_stars(:,2);     %hodnoty jasnosti = y-ova os
X = [ones(size(x),1),log10(x)];     %matica derivacii modelovej funkcie podla parametrov
W = eye(size(x,1));     %jednotkova diagonalna matica reprezentujuca vahy
V = X'*W*X;
U = X'*W*y;
b = inv(V)*U;
z = b(1)+b(2)*log10(x);     %fit
plot(x,y,'*',x,z,'r')
xlabel('poradie')
ylabel('hviezdna velkost [mag]')
legend('data','fit','Location','SouthEast')

%uloha 2
[HJD,M4,M0,M1,M2,M3,E4,E0,E1,E2,E3,G,F] = textread('V2509_Sgr.dat','%f %f %f %f %f %f %f %f %f %f %f %s %f');
fi = HJD/1.087 - floor(HJD/1.087);  %fazova funkcia
fi = [fi-1;fi;fi+1];    %opakovanie fazoveho intervalu
E2 = [E2;E2;E2];
M2 = [M2;M2;M2];
figure
for i = 1:6
    plot(fi,M2,'*','MarkerSize',2)
    n = size(M2,1);
    X = [ones(n,1),cos(2*pi*fi),sin(2*pi*fi),cos(4*pi*fi),sin(4*pi*fi),cos(6*pi*fi),sin(6*pi*fi)];
    W = diag(E2.^(-2)/mean(E2.^(-2)));
    V = X'*W*X;
    U = X'*W*M2;
    b = inv(V)*U;
    t = -1:0.01:2;
    u = b(1)+b(2)*cos(2*pi*t)+b(3)*sin(2*pi*t)+b(4)*cos(4*pi*t)+b(5)*sin(4*pi*t)+b(6)*cos(6*pi*t)+b(7)*sin(6*pi*t);
    hold on
    plot(t,u,'r')
    axis ij
    xlabel('faza')
    ylabel('hviezdna velkost [mag]')
    hold off
    sigma = sqrt((M2'*W*M2-b'*U)/(n-size(X,2)));    %standardna odchylka
    error = find(abs(M2-X*b)>3.5*sigma);    %odlahle body
    fi(error) = [];     %odstranenie odlahlych bodov
    M2(error) = [];
    E2(error) = [];
    pause       %stlacenim klavesy v prikazovom riadku sa spusti dalsi cyklus
end