load asu.txt   %asu.txt = stiahnuty datovy subor
t = asu(:,2);   %HJD
y = asu(:,3);   %ETV
w = asu(:,4).^(-2);   %vahy
plot(t,y,'.')
M0 = 55950;   %hodnota nuloveho HJD; pre pouzivanu fitovaciu funkciu je to hodnota vrcholu (idealne v strednej oblasti HJD)
xlabel('HJD [d]')
ylabel('intensity [mag]')
figure
Sigma = [];   Period = [];
for P = 60:0.1:80   %hodnotu priblizne odhadujeme z predchadzajuceho grafu, mozme menit vzorkovaciu frekvenciu
    phi = (t-M0)/P-floor((t-M0)/P);   %fazova funkcia
    b = [0,50,0]';   %pociatocne hodnoty parametrov
    for i = 1:5   %pocet cyklov nelinearnej regresie
        f0 = b(1);   A = b(2);   phi0 = b(3);
        X = [phi./phi,cos(2*pi*(phi-phi0)),A*2*pi*sin(2*pi*(phi-phi0))];
        f = f0+A*cos(2*pi*(phi-phi0));   %funkcne hodnoty fitu
        Dy = y-f;
        V = X'*(repmat(w,1,3).*X);   %pocetne zjednoduseny vyraz v porovnani s X'*W*X
        U = X'*(w.*Dy);
        Db = inv(V)*U;
        b = Db+b;
    end
    sigma = sqrt(Dy'*(w.*Dy)/(size(X,1)-3)/mean(w));   %standardna odchylka
    if sigma < min(Sigma)   %vykresli graf s najlepsou periodou
        plot(phi,y,'.',phi,f,'r.')
        xlabel('phase')
        ylabel('intensity [mag]')
        title(['P = ', num2str(P), ' d'])
    end
    Sigma = [Sigma;sigma];
    Period = [Period;P];
end
figure
plot(Period,Sigma)
xlabel('period [d]')
ylabel('standard deviation')