clc;
%Pontos visados
%Clérigos
clr_M = -40420.448;
clr_P = 164167.927;
clr_H = 141.13;
%Direta Progressiva (dp)
clr_ah_dp = 102.8444;
clr_av_dp = 97.7522;
%Inversa Regressiva (ir)
clr_ah_ir = 302.8566;
clr_av_ir = 302.2442;
%Inversa Progressiva (ip)
clr_ah_ip = 302.8561;
clr_av_ip = 302.2467;
%Direta Regressiva (dr)
clr_ah_dr = 102.8617;
clr_av_dr = 97.7503;
%Câmara municipal do porto
cm_M = -40088.726;
cm_P = 164636.884;
cm_H = 152.33;
%Direta Progressiva (dp)
cm_ah_dp = 129.6582;
cm_av_dp = 97.4767;
%Inversa Regressiva (ir)
cm_ah_ir = 329.6599;
cm_av_ir = 302.5056;
%Inversa Progressiva (ip)
cm_ah_ip =329.6628;
cm_av_ip =302.3095;
%Direta Regressiva (dr)
cm_ah_dr = 129.6708;
cm_av_dr = 97.4755;
%igreja da lapa
il_M = -40209.736;
il_P = 165426.432;
il_H = 176.03;
%Direta Progressiva (dp)
il_ah_dp = 137.6340 ;
il_av_dp = 97.6790;
%Inversa Regressiva (ir)
il_ah_ir = 337.6339;
il_av_ir = 302.3148;
%Inversa Progressiva (ip)
il_ah_ip = 337.6455;
il_av_ip = 302.3175;
%Direta Regressiva (dr)
il_ah_dr = 137.6501;
il_av_dr = 97.6786;
%Jornal de Noticias
jn_M = -39934.342;
jn_P = 165011.012;
jn_H = 170.79;
%Direta Progressiva (dp)
jn_ah_dp = 141.7560;
jn_av_dp = 97.2225;
%Inversa Regressiva (ir)
jn_ah_ir = 341.7596;
jn_av_ir = 302.7853;
%Inversa Progressiva (ip)
jn_ah_ip = 341.7649;
jn_av_ip = 302.7908;
%Direta Regressiva (dr)
jn_ah_dr = 141.7680;
jn_av_dr = 97.2106;
%Depósito de água
da_M = -39336.276;
da_P = 165413.16;
da_H = 185.48;
%Direta Progressiva (dp)
da_ah_dp = 163.9711;
da_av_dp = 97.1290;
%Inversa Regressiva (ir)
da_ah_ir = 363.9730;
da_av_ir = 302.8053;
%Inversa Progressiva (ip)
da_ah_ip = 363.9772;
da_av_ip = 302.8023;
%Direta Regressiva (dr)
da_ah_dr = 163.9848;
da_av_dr = 97.1911;
%Igreja do Bonfim
ib_M = -38614.9280;
ib_P = 164792.1040;
ib_H = 162.20;
%Direta Progressiva (dp)
ib_ah_dp = 194.7617;
ib_av_dp = 97.3887;
%Inversa Regressiva (ir)
ib_ah_ir = 394.7702;
ib_av_ir = 302.5992;
%Inversa Progressiva (ip)
ib_ah_ip = 394.7724;
ib_av_ip = 302.5994;
%Direta Regressiva (dr)
ib_ah_dr = 194.7772;
ib_av_dr = 97.3917;
%Gondomar(v.g.)
g_M = -33142.6520;
g_P = 163536.09;
g_H = 203.27;
%Direta Progressiva (dp)
g_ah_dp = 260.0913;
g_av_dp = 98.9460;
%Inversa Regressiva (ir)
g_ah_ir = 60.0908;
g_av_ir = 301.0498;
%Inversa Progressiva (ip)
g_ah_ip = 60.0992;
g_av_ip = 301.0477;
%Direta Regressiva (dr)
g_ah_dr = 260.0978;
g_av_dr = 98.9363;
%
%3 Pontos- Clérigos, Igreja do bonfim e jornal de noticias-
dir_clr = ((clr_ah_dp+clr_ah_ir-200+clr_ah_ip-200+clr_ah_dr)/4)*(180/200);
dir_cm = ((cm_ah_dp+cm_ah_ir-200+cm_ah_ip-200+cm_ah_dr)/4)*(180/200);
dir_il = ((il_ah_dp+il_ah_ir-200+il_ah_ip-200+il_ah_dr)/4)*(180/200);
dir_jn = ((jn_ah_dp+jn_ah_ir-200+jn_ah_ip-200+jn_ah_dr)/4)*(180/200);
dir_da = ((da_ah_dp+da_ah_ir-200+da_ah_ip-200+da_ah_dr)/4)*(180/200);
dir_ib = ((ib_ah_dp+ib_ah_ir-200+ib_ah_ip-200+ib_ah_dr)/4)*(180/200);
dir_g = ((g_ah_dp+g_ah_ir-200+g_ah_ip-200+g_ah_dr+400)/4)*(180/200);
%Calculo rumos BA e BC (ponto A - Câmara do Porto, B- Depósito de água e C-
%Igreja do bonfim)
R_BA = atand((clr_M-jn_M)/(clr_P-jn_P))+180;
R_BC = atand((ib_M-jn_M)/(ib_P-jn_P))+180;
%Calculo do angulo phi
a_phi=R_BA-R_BC;
%Calculo distâncias horizontais
d_BA = sqrt((clr_M-jn_M)^2+(clr_P-jn_P)^2);
d_BC = sqrt((ib_M-jn_M)^2+(ib_P-jn_P)^2);
%Calculo do R
a_alpha = dir_jn - dir_cm;
a_beta = dir_ib - dir_jn;
R = 360-(a_alpha+a_beta+a_phi);
%Calculo do S
S=(d_BC+sind(a_alpha))/(d_BA*sind(a_beta));
%Calculo do gama
a_gama = atand(sind(R)/(S+cosd(R)))+90;
%Calculo tetha
a_theta=R-a_gama;
%calculo phi1
a_phi1= 180-(a_theta+a_alpha);
%Calculo distancia Camara do Porto(A) ao ponto P(local observação)
d_AP=d_BA*((sind(a_phi1))/sind(a_alpha));
%Calculo Rumo AP
RAP = R_BA-180+a_theta;
%Coordenadas aproximadas M0 e P0
M0 = clr_M + d_AP*sind(RAP);
P0 = clr_P + d_AP*cosd(RAP);
%%
h_medido = [clr_H, jn_H, ib_H];
%Distancia ao ponto aproximado
d_clr = sqrt((M0-clr_M)^2+(P0-clr_P)^2);
d_cm = sqrt((M0-cm_M)^2+(P0-cm_P)^2);
d_il = sqrt((M0-il_M)^2+(P0-il_P)^2);
d_jn = sqrt((M0-jn_M)^2+(P0-jn_P)^2);
d_da = sqrt((M0-da_M)^2+(P0-da_P)^2);
d_ib = sqrt((M0-ib_M)^2+(P0-ib_P)^2);
d_g = sqrt((M0-g_M)^2+(P0-g_P)^2);
%Rumos ao ponto aproximado
r_clr = atand((clr_M-M0)/(clr_P-P0));
r_cm = atand((cm_M-M0)/(cm_P-P0));
r_il = atand((il_M-M0)/(il_P-P0))+360;
r_jn = atand((jn_M-M0)/(jn_P-P0));
r_da = atand((da_M-M0)/(da_P-P0));
r_ib = atand((ib_M-M0)/(ib_P-P0))+360;
r_g = atand((g_M-M0)/(g_P-P0))+360;
%Matriz dos coeficientes
A = [-(clr_P-P0)/d_clr^2*180/pi (clr_M-M0)/d_clr^2*180/pi -1;
-(cm_P-P0)/d_cm^2*180/pi (cm_M-M0)/d_cm^2*180/pi -1;
-(il_P-P0)/d_il^2*180/pi (il_M-M0)/d_il^2*180/pi -1;
-(jn_P-P0)/d_jn^2*180/pi (jn_M-M0)/d_jn^2*180/pi -1;
-(da_P-P0)/d_da^2*180/pi (da_M-M0)/d_da^2*180/pi -1;
-(ib_P-P0)/d_ib^2*180/pi (ib_M-M0)/d_ib^2*180/pi -1;
-(g_P-P0)/d_g^2*180/pi (g_M-M0)/d_g^2*180/pi -1];
%Vetor b
b = [dir_clr-r_clr;
dir_cm-r_cm;
dir_il-r_il;
dir_jn-r_jn;
dir_da-r_da;
dir_ib-r_ib;
dir_g-r_g];
%Matriz dos pesos
W=eye(7);
%Matriz equações normais
N=A'*W*A ;
%Calculo correções coordenadas
delta = inv(N) * A' * W * b;
%coordenadas compensadas
M0_f=M0+delta(1);
P0_f=P0+delta(2);
disp(['Coordenadas finais M0_f = ', num2str(M0_f)]);
disp(['Coordenadas finais P0_f = ', num2str(P0_f)]);
%%
% Cálculo da cota corrigida
R_T = 6371000; % Raio da Terra em metros
k = 0.14; % Coeficiente de refração atmosférica
% Distâncias horizontais correspondentes
d_P0_N1 = sqrt((M0_f - clr_M)^2 + (P0_f - clr_P)^2);
d_P0_N2 = sqrt((M0_f - jn_M)^2 + (P0_f - jn_P)^2);
d_P0_N3 = sqrt((M0_f - ib_M)^2 + (P0_f - ib_P)^2);
dist = [d_P0_N1, d_P0_N2, d_P0_N3];
% Correções de curvatura e refração
cc = @(dist) (dist^2) / (2 * R_T);
rc = @(C) k * C;
% Cálculo das cotas corrigidas
h_corrigido = zeros(1, 3);
for i = 1:length(h_medido)
C = cc(dist(i));
R = rc(C);
h_corrigido(i) = h_medido(i) - (C + R);
end
% Média das cotas corrigidas (valor mais provável)
h_mean = mean(h_corrigido);
% Resultados das cotas corrigidas
% Exibir resultados intermediários
fprintf('Correções: ΔM = %.4f, ΔP = %.4f, R0 = %.4f\n', delta(1), delta(2), delta(3));
fprintf('Coordenadas compensadas: M0_f = %.4f, P0_f = %.4f\n', M0_f, P0_f);
% Calcular a matriz de covariância e precisão da solução
residuo = b - A * delta;
sigma2 = sum(residuo.^2) / (length(b) - length(delta)); % Variância residual média
cov_matrix = sigma2 * inv(A' * A);
% Precisão das coordenadas compensadas (desvios padrão)
std_devs = sqrt(diag(cov_matrix));
% Exibir matriz de covariância e desvios padrão
disp('Matriz de covariância dos parâmetros:');
disp(cov_matrix);
disp('Precisão das coordenadas compensadas (desvios padrão):');
disp(std_devs);
this is the progress but i still feel like its not good enough or its wrong even.
