load SQ.dat
MJD = SQ(:,1);   %modifikovane julianske datum
int = SQ(:,2);   %intenzita vo filtri V
plot(MJD,int,'.')
t0 = [0];   %pociatocna hodnota odhadu polohy centra zakrytu
for i = 1:size(MJD,1)
    if int(i)>8.4 & abs(MJD(i)-t0(end))>10
        t0 = [t0;MJD(i)];   %odhadovane polohy centier vsetkych zakrytov
    end
end
t0(1) = [];   %umelo dodanu prvu hodnotu musime odstranit
f0 = 8.15;   A = 0.3;   s = 0.7;
b = [f0;A;s;t0];   %pociatocne parametre modelovej funkcie
t = min(MJD):0.1:max(MJD);   %spojity casovy vektor pre vykreslenie fitovanej funkcie ft
for j = 1:10
    f0 = b(1);   A = b(2);   s = b(3);   t0 = b(4:end);
    X1 = ones(size(MJD));   %derivacia podla f0
    X2 = 0;   X3 = 0;   X4 = [];   f = f0;   ft = f0;
    for i = 1:size(t0,1)
        X2 = X2+exp(-(MJD-t0(i)).^2/2/s^2);   %derivacia podla A
        X3 = X3+A*exp(-(MJD-t0(i)).^2/2/s^2).*(MJD-t0(i)).^2/s^3;   %derivacia podla s
        X4 = [X4,A*exp(-(MJD-t0(i)).^2/2/s^2).*(MJD-t0(i))/s^2];   %derivacie podla t0i
        f = f+A*exp(-(MJD-t0(i)).^2/2/s^2);   %funkcna hodnota v casoch MJD
        ft = ft+A*exp(-(t-t0(i)).^2/2/s^2);   %funkcna hodnota v casoch t
    end
    plot(MJD,int,'.',t,ft,'r')
    axis ij
    xlabel('MJD [d]')
    ylabel('intensity [mag]')
    title('SQ Tau lightcurve')
    X = [X1,X2,X3,X4];
    DY = int-f;
    V = X'*X;
    U = X'*DY;
    Db = inv(V)*U;   %korekcny clen Delta b: b_j = b_{j-1} + Delta b
    b = b+Db;
    sigma = sqrt(DY'*DY/(size(X,1)-size(X,2)));
    db = sqrt(sigma^2*diag(inv(V)));   %chyba urcenia parametra b
end