h = [2.370; 4.718; -8.224; 2.909; -1.810];
m = [4; 4; 2; 2; 4]*1e-3;
H = [100.127; 115.430];

A = [
   1  0  0
  -1  1  0
   0  1  0
   0  1 -1
   1  0 -1
  ] ;

P = diag(m.^-2) ;

L = [ -H(1) ; 0 ; -H(2) ; 0 ; 0 ] - h ;

X = -inv(A'*P*A)*(A'*P*L);

printf("\nWyrównane wysokości:\n")
printf("%10.4f\n", X)

V = A*X + L ;
hw = h + V;

printf("\nWyrównane przewyższenia:\n")
printf("%10.4f\n", hw)


[n,r]  = size(A);
sigma0 = sqrt(V'*P*V/(n-r));

Cx = sigma0**2*(A'*P*A)^(-1);
mH = sqrt(diag(Cx))

rho = [];
for i = 1:size(Cx,1)-1
  for j = i+1:size(Cx,1)
    rho = [rho ; Cx(i,j)/(mH(i)*mH(j))];
  end
end
rho

Cv = sigma0^2 * (inv(P) - A*inv(A'*P*A)*A');
mv = sqrt(diag(Cv))

Ch = sigma0**2 * (A*inv(A'*P*A)*A');
mh = sqrt(diag(Ch))

D = [-1 1 0] ;
Cu = D * Cx * D';
sqrt(Cu)

D = [1 1 1 0 0 ; 1 -1 0 0 0 ];
Cu= D*Ch*D';
mu = sqrt(diag(Cu))
