Q = load('cont0109a.txt');
n = size(Q,1);
x = Q(:,1);     %vlnova dlzka
y = Q(:,2);     %intenzita
b0 = [1,-0.7,4340,3]';
for i = 1:5
    G0 = b0(1);     %intenzita kontinua
    A = b0(2);      %amplituda
    x0 = b0(3);     %stred ciary
    sigma = b0(4);  %polosirka ciary
    G = G0+A*exp(-(x-x0).^2/2/sigma^2);     %prekladana Gaussova funkcia
    X = [ones(n,1),exp(-(x-x0).^2/2/sigma^2),A*exp(-(x-x0).^2/2/sigma^2).*(x-x0)/sigma^2,A*exp(-(x-x0).^2/2/sigma^2).*(x-x0).^2/sigma^3];
    dY=y-G;
    W=eye(n);
    V=X'*W*X;
    U=X'*W*dY;
    db=inv(V)*U;
    b0 = db+b0;
    s=sqrt((y'*W*y-db'*U)/(n-4));
    plot(x,y,'b',x,G,'r')
    pause   %pri zobrazeni pokracujeme enterom
end
s

b0 = [1,-0.7,4340,3]';
for i = 1:5
    L0 = b0(1);
    A = b0(2);
    x0 = b0(3);
    gamma = b0(4);
    L = L0+A*gamma./((x-x0).^2+gamma^2);    %prekladana Lorentzova funkcia
    X = [ones(n,1),gamma./((x-x0).^2+gamma^2),A*2*gamma*(x-x0)./((x-x0).^2+gamma^2).^2,A*((x-x0).^2-gamma^2)./((x-x0).^2+gamma^2).^2];
    dY=y-L;
    W=eye(n);
    V=X'*W*X;
    U=X'*W*dY;
    db=inv(V)*U;
    b0 = db+b0;
    s=sqrt((y'*W*y-db'*U)/(n-4));
    plot(x,y,'b',x,G,'r',x,L,'g')
    legend('data', 'Gauss', 'Lorentz', 'Location', 'SouthEast')
    xlabel('vlnova dlzka [A]')
    ylabel('relativna intenzita')
    pause
end
s