clc;
clear all;

format short

disp (' ***********************************************************')
disp (' Ajustamento de uma linha fechada de nivelamento geometrico ')
disp (' ***********************************************************')

% as equacoes (lineares, no caso de nivelamento geometrico)que relacionam 
% os parametros com as observacoes sao do tipo AX+V+W=0, onde A e a matriz 
% dos coeficientes dos parametos, X e o vector dos parametros, V e o vector
% dos resíduos e W e o vector dos termos independentes (cotas e desniveis)

% vector com o nome dos pontos

pontos=cellstr(['MN1';'  A';'  B';' B1';'  C';'  D';'  E';'  F';'MN2']);

% vector dos desniveis

l=[0.836;-0.267;-0.699;0.544;0.701;0.637;-0.345;0.568];

l

pause

somatorio_desniveis=l(1,1)+l(2,1)+l(3,1)+l(4,1)+l(5,1)+l(6,1)+l(7,1)+l(8,1);

erro_de_fecho=196.251-198.219+somatorio_desniveis

pause

% matriz A dos coeficientes dos parametros

A=zeros(8,7);

A(1,1)=-1;
A(2,1)=1;
A(2,2)=-1;
A(3,2)=1;
A(3,3)=-1;
A(4,3)=1;
A(4,4)=-1;
A(5,4)=1;
A(5,5)=-1;
A(6,5)=1;
A(6,6)=-1;
A(7,6)=1;
A(7,7)=-1;
A(8,7)=1;


A

pause

% vector W dos termos independentes (desniveis e/ou cotas dos pontos 
% de partida e chegada), m

W(1,1)=196.251+0.836;
W(2,1)=-0.267;
W(3,1)=-0.699;
W(4,1)=0.544;
W(5,1)=0.701;
W(6,1)=0.637;
W(7,1)=-0.345;
W(8,1)=0.568-198.219;

W

pause

% vector das distancias entre miras, m

dist(1,1)=66.9;
dist(2,1)=59.1;
dist(3,1)=35.7;
dist(4,1)=45.2;
dist(5,1)=52.8;
dist(6,1)=65.3;
dist(7,1)=26.1;
dist(8,1)=15.0;

dist

pause

% precisão corrente (slide 84)
tolerancia=2.6*4.5/206265*sqrt(dist(1,1)^2+dist(2,1)^2+dist(3,1)^2+dist(4,1)^2+dist(5,1)^2+dist(6,1)^2+dist(7,1)^2+dist(8,1)^2)

% pag 112 Topografia, conceitos e aplicacoes, precisão corrente
tolerancia1=12*sqrt((dist(1,1)+dist(2,1)+dist(3,1)+dist(4,1)+dist(5,1)+dist(6,1)+dist(7,1)+dist(8,1))/1000)/1000

pause

% desvio padrao do Wild NA2: 0.7 mm/km = 0.0007 m/km

sigma_NA2=0.0007;
%sigma_NA2=0.007;

sigma=zeros(8,8);

sigma(1,1)=sigma_NA2*sqrt(2*dist(1,1)/1000);
sigma(2,2)=sigma_NA2*sqrt(2*dist(2,1)/1000);
sigma(3,3)=sigma_NA2*sqrt(2*dist(3,1)/1000);
sigma(4,4)=sigma_NA2*sqrt(2*dist(4,1)/1000);
sigma(5,5)=sigma_NA2*sqrt(2*dist(5,1)/1000);
sigma(6,6)=sigma_NA2*sqrt(2*dist(6,1)/1000);
sigma(7,7)=sigma_NA2*sqrt(2*dist(7,1)/1000);
sigma(8,8)=sigma_NA2*sqrt(2*dist(8,1)/1000);

% matriz sigma2 das variancias-covariancias das observacoes (desniveis), m^2
% (observacoes nao correlacionadas)

sigma2=sigma^2;

% variancia a priori por unidade de peso

sigma02=1;

% matriz cofactor QL das observacoes (desniveis): matriz sigma2 escalada por sigma02

QL=sigma2/sigma02;

% matriz dos pesos dos desniveis

P=inv(QL);

P

pause

% calculo do ajustamento

At=transp(A);

N=At*P*A;

U=At*P*W;

% parametros ajustados

X=-inv(N)*U;

cotas=X;

cotas

pause

% matriz cofactor dos parametros

Qx=inv(N);

Qx;

tracoQx=Qx(1,1)+Qx(2,2)+Qx(3,3)+Qx(4,4)+Qx(5,5)+Qx(6,6)+Qx(7,7)

% precisao das cotas

for i=1:7
    precisao_cotas(i,1)=sqrt(Qx(i,i));
end

precisao_cotas

pause

% residuos das observacoes

v=-A*X-W;

v

% matriz cofactor dos residuos

Qv=(inv(P)-A*inv(N)*At);

Qv

pause

% observacoes ajustadas

La(1,1)=l(1,1)+v(1,1);
La(2,1)=l(2,1)+v(2,1);
La(3,1)=l(3,1)+v(3,1);
La(4,1)=l(4,1)+v(4,1);
La(5,1)=l(5,1)+v(5,1);
La(6,1)=l(6,1)+v(6,1);
La(7,1)=l(7,1)+v(7,1);
La(8,1)=l(8,1)+v(8,1);

La

erro_fecho1=196.251-198.219+La(1,1)+La(2,1)+La(3,1)+La(4,1)+La(5,1)+La(6,1)+La(7,1)+La(8,1);

erro_fecho1


% matriz cofactor das observacoes ajustadas

QLa=inv(P)-Qv;

% variancia a posteriori por unidade de peso

s02=transp(v)*P*v/(8-7)

% distribuicao do qui-quadrado a 0.05%

qui(1,1)=3.841;
qui(2,1)=5.991;
qui(3,1)=7.815;
qui(4,1)=9.488;
qui(5,1)=11.071;

% teste da razao de variancias

y=s02*(13-12)/sigma02;

if (y < qui(8-7,1)) 
% sigma02 e s02 estatisticamente semelhantes    
    'Aceitar o ajustamento.'
else
    'Rejeitar o ajustamento.'
end


































