clc
clear all
close all

format long

disp('  ')
disp('Exemplo slide 89 IM1')
disp('  ')

% valores observados

ab=136.45; %em m
sigma_ab=0.018; %em m

alfa=(44+12/60+34/3600)*pi/180; %em radiano
sigma_alfa=8.6/3600*pi/180; %em radiano

beta=(39+26/60.0+56/3600)*pi/180; %em radiano
sigma_beta=11.3/3600*pi/180; %em radiano

zA=(81+47/60+13/3600)*pi/180; %em radiano
sigma_zA=4.1/3600*pi/180; %em radiano

zB=(84+9.0/60+50/3600)*pi/180; %em radiano
sigma_zB=5.1/3600*pi/180; %em radiano

HA=1298.65; %em m
sigma_HA=0.006; %em m

HB=1301.53; %em m
sigma_HB=0.004; %em m

aiA=5.25; %em m
sigma_aiA=0.005; %em m

aiB=5.18; %em m
sigma_aiB=0.005; %em m

% resolucao do triangulo pela analogia dos senos

ac=ab*sin(beta)/(sin(alfa+beta));
bc=ab*sin(alfa)/(sin(alfa+beta));

% calculo da altitude do ponto C

HCA=HA+aiA+ac/tan(zA);
HCB=HB+aiB+bc/tan(zB);

HC=(HCA+HCB)/2;

fprintf('altitude do ponto C = %8.3f m\n', HC)

syms HC1 HA1 HB1 aiA1 aiB1 alfa1 beta1 zA1 zB1 ab1
HC1=(HA1+aiA1+ab1*sin(beta1)/(tan(zA1)*sin(alfa1+beta1))+HB1+aiB1+ab1*sin(alfa1)/(tan(zB1)+sin(alfa1+beta1)))/2;
% calcular as derivadas parciais da funcao HC1 em relacao as observacoes
d_HC1_HA1=diff(HC1,HA1);
d_HC1_HB1=diff(HC1,HB1);
d_HC1_aiA1=diff(HC1,aiA1);
d_HC1_aiB1=diff(HC1,aiB1);
d_HC1_alfa1=diff(HC1,alfa1);
d_HC1_beta1=diff(HC1,beta1);
d_HC1_zA1=diff(HC1,zA1);
d_HC1_zB1=diff(HC1,zB1);
d_HC1_ab1=diff(HC1,ab1);

J(1,1)=double(subs(d_HC1_HA1,[HA1 HB1 aiA1 aiB1 alfa1 beta1 zA1 zB1 ab1],[HA HB aiA aiB alfa beta zA zB ab]));
J(1,2)=double(subs(d_HC1_HB1,[HA1 HB1 aiA1 aiB1 alfa1 beta1 zA1 zB1 ab1],[HA HB aiA aiB alfa beta zA zB ab]));
J(1,3)=double(subs(d_HC1_aiA1,[HA1 HB1 aiA1 aiB1 alfa1 beta1 zA1 zB1 ab1],[HA HB aiA aiB alfa beta zA zB ab]));
J(1,4)=double(subs(d_HC1_aiB1,[HA1 HB1 aiA1 aiB1 alfa1 beta1 zA1 zB1 ab1],[HA HB aiA aiB alfa beta zA zB ab]));
J(1,5)=double(subs(d_HC1_alfa1,[HA1 HB1 aiA1 aiB1 alfa1 beta1 zA1 zB1 ab1],[HA HB aiA aiB alfa beta zA zB ab]));
J(1,6)=double(subs(d_HC1_beta1,[HA1 HB1 aiA1 aiB1 alfa1 beta1 zA1 zB1 ab1],[HA HB aiA aiB alfa beta zA zB ab]));
J(1,7)=double(subs(d_HC1_zA1,[HA1 HB1 aiA1 aiB1 alfa1 beta1 zA1 zB1 ab1],[HA HB aiA aiB alfa beta zA zB ab]));
J(1,8)=double(subs(d_HC1_zB1,[HA1 HB1 aiA1 aiB1 alfa1 beta1 zA1 zB1 ab1],[HA HB aiA aiB alfa beta zA zB ab]));
J(1,9)=double(subs(d_HC1_ab1,[HA1 HB1 aiA1 aiB1 alfa1 beta1 zA1 zB1 ab1],[HA HB aiA aiB alfa beta zA zB ab]));

%disp ('matriz de variancias-covariancias das observacoes')
Sigma_obs=eye(9,9);
Sigma_obs(1,1)=sigma_HA^2;
Sigma_obs(2,2)=sigma_HB^2;
Sigma_obs(3,3)=sigma_aiA^2;
Sigma_obs(4,4)=sigma_aiB^2;
Sigma_obs(5,5)=sigma_alfa^2;
Sigma_obs(6,6)=sigma_beta^2;
Sigma_obs(7,7)=sigma_zA^2;
Sigma_obs(8,8)=sigma_zB^2;
Sigma_obs(9,9)=sigma_ab^2;


Sigma_HC=J*Sigma_obs*J';

incerteza_HC=sqrt(Sigma_HC);

fprintf('incerteza na altitude do ponto C = %6.3f m\n', incerteza_HC)

disp('  ')

