HR = 1;
h = [1.04 ; 2.05 ; 3.03];
m = [ 2 ; 2 ; 4]/1e2;

%macierz współczynników przy niewiadomych
A = [ 1 0 ; -1 1; 0 1;];
% macierz wag
P = diag(m.^(-2));
% wektor wyrazów wolnych
L = [-HR ; 0 ; -HR ] - h;

% policzenie wysokości wyrównanych
X = -inv(A'*P*A)*(A'*P*L);
% policzenie poprawek
V = A*X+L;
% policzenie obserwacji wyrównanych
hw = h + V;

% kontrola => wektor 2x1 zer
A'*P*V;

% błąd średni typowego spostrzeżenia
sigma = sqrt(V'*P*V/(size(A,1)-size(A,2)));

% macierz kowariancyjna wektora niewiadomych
Cx = sigma^2*inv(A'*P*A);
% błędy wysokości
mH = sqrt(diag(Cx));

korelacje_H1H2 = Cx(1,2)/(mH(1)*mH(2));% musi być <-1,1>;

% macierz kowariancyjna wektora poprawek
CV = sigma^2*(inv(P) - A*inv(A'*P*A)*A');
% błędy poprawek
mV = sqrt(diag(CV));

% macierz kowariancyjna wektora obserwacji
Ch = sigma^2*(A*inv(A'*P*A)*A');
% błędy obserwacji po wyrównaniu
mh = sqrt(diag(Ch));

% błąd funkcji u = H2-H1
D = [ -1 1 ];
Cu = D*Cx*D';
mu = sqrt(Cu);

%ostatni punkt ćwiczenia
D = [ 0  2 ; 2 -1];% lub +1
Cu = D*Cx*D';
mu = sqrt(diag(Cu))
% lub
D = [ 1 1 1 ; 1 -1 0];
Cu = D*Ch*D';
mu = sqrt(diag(Cu))

