clc
close all
clear all

sigma02=1;
nsign=0.05;

% desniveis m
% epoca 1
% l(7) passa de -0.378 para -0.280
l=[2.474;-1.250;-1.241;-2.233;3.117;-2.115;-0.378;-3.094;2.822];
%l=[2.474;-1.250;-1.241;-2.233;3.117;-2.115;-0.280;-3.094;2.822];

% epoca 2
%l=[2.470;-1.253;-1.240;-2.235;3.111;-2.114;-0.278;-3.094;2.819];

% comprimento desniveis km
km=[4.1;2.2;2.4;6.0;3.6;4.0;2.2;2.3;3.5];

% incerteza desniveis m
for i=1:9
    sigma_l(i,i)=sqrt(km(i))/100;
    i=i+1;
end

sigma_l=sigma_l';

Q_l=sigma_l.^2/sigma02;

W_l=inv(Q_l)

n=9;
n0=5;
r=n-n0;

A=zeros(4,9);
f=zeros(4,1);

% metodo condicional: 4 eqs que relacionam os 9 desniveis, Av=f
% I:   l1+v1+l2+v2+l3+v3=0       => v1+v2+v3=-l1-l2+l3
% II:  l2+v2-l6-v6-l5-v5-l4-v4=0 => v2-v4-v5-v6=-l2+l4+l5+l6
% III: l3+v3+l6+v6-l7-v7-l8-v8=0 => v3+v6-v7-v8=-l3-l6+l7+l8
% IV:  l5+v5+l7+v7-l9-v9=0       => v5+v7-v9=-l5-l7+l9

A=[1 1 1 0 0 0 0 0 0;0 1 0 -1 -1 -1 0 0 0;0 0 1 0 0 1 -1 -1 0;0 0 0 0 1 0 1 0 -1]
f=[-l(1)-l(2)-l(3);-l(2)+l(4)+l(5)+l(6);-l(3)-l(6)+l(7)+l(8);-l(5)-l(7)+l(9)]

k=inv(A*inv(W_l)*A')*f;

v=inv(W_l)*A'*k;


Q_v=Q_l*A'*inv(A*inv(W_l)*A')*A*Q_l;

l_ajust=l+v;

Q_l_ajust=Q_l-Q_v;

fprintf('Residuos das observações (m): \n');
for i=1:9
    fprintf('      %+.5f\n', v(i));
    i=i+1;
end

disp ('   ')

s02=v'*W_l*v/r

disp 'Teste da razao de variancias:'
disp '                             '

quiM=icdf('chi2',1-nsign/2,r);
quim=icdf('chi2',nsign/2,r);
        
if (quim < s02*r/sigma02 & s02*r/sigma02 < quiM)
    disp '  --> Aceitar o ajustamento'
    op2=0;
else
    disp '  --> Rejeitar o ajustamento'  
    op2=1;
end

if (op2==0)
    sigma_l_ajust=s02*Q_l_ajust   
    incerteza_global=trace(sqrt(sigma_l_ajust))
end    
               
disp('    ')
disp('Fim do programa.')     
