[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];   %chyba (error)
M2 = [M2;M2;M2];   %jasnost (magnitude)
%vyberame kombinaciu jasnost + chyba, alternativne mozme zvolit M0+E0, M1+E1 alebo M4+E4
matrix = sortrows([fi,E2,M2],1);
fi = matrix(:,1);   E2 = matrix(:,2);   M2 = matrix(:,3);
%zoradujeme vektory fi, E2, M2 podla fazy, aby v prikaze plot spajali linie fazovo susedne body
for i = 1:5
    plot(fi,M2,'*')
    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)];
    %stupen harmonickeho polynomu urcime na zaklade vlastneho uvazenia
    W = diag(E2.^(-2));
    V = X'*W*X;
    U = X'*W*M2;
    b = inv(V)*U;
    yp = X*b;
    hold on
    plot(fi,yp,'r')
    axis ij
    xlabel('faza')
    ylabel('hviezdna velkost [mag]')
    chi2 = M2'*W*M2-b'*U   %chi kvadrat
    g = size(X,2);
    chi2_mu = chi2/(n-g)   %modifikovany chi kvadrat
    sigma = sqrt(chi2_mu/mean(diag(W)))
    dyp = sqrt(chi2_mu*diag(X*inv(V)*X'));
    db = sqrt(chi2_mu*diag(inv(V)));
    plot(fi,yp+3*dyp,'r:',fi,yp-3*dyp,'r:')   %interval spolahlivosti urcenia fitu s pravdepodobnostou 3*sigma (99.7%)
    hold off
    error = find(abs(M2-X*b)>3.5*sigma);   %odlahle body
    fi(error) = [];   %odstranenie odlahlych bodov
    M2(error) = [];
    E2(error) = [];
    pause
end