Vyrovnání oboustraně připojeného a orientovaného polygonového pořadu pomocí MNČ
zpracoval: Zdeněk Nejedlý
dne: 27.12.2002
Ovládání programu | Příklad | Zdrojový kód
Soubor "plg.m" je třeba umístit do aktuálního pracovního adresáře (Current Directory) nebo do některého z adresářů uvedených v [File / Set Path...] jinak ho Matlab nenajde a nebude s ním počítat. Do stejného adresáře je vhodné umístit i seznam souřadnic a zápisník s měřenými daty (tedy pokud nedáš přednost ručnímu zadání hodnot)
Ovládání programu
Program se spustí příkazem plg('seznam souřadnic','zápisník měřených dat'), konkrétně třeba plg('ss.txt','zap.txt'). Jméno soubodu je třeba zadat i s apostrofy ('). Pokud jsou vstupní soubory v jiném adresáři než program, je třeba zadat jejich úplnou cestu, například plg('g:/ss.txt','c:\data/zap.txt')
Pokud se program spustí bez parametrů, objeví se nabídka s volbou vstupu
Nacist data ze souboru - 1
Zadat data rucne - 2
Typ vstupu :
Po stisknutí 1 můžeš zadat vstupní soubory, ...
Nacist data ze souboru - 1
Zadat data rucne - 2
Typ vstupu : 1
Nacteni dat ze souboru
(jmeno souboru zadavej ve tvaru jmeno.pripona (bez'))
seznam souradnic = ss.txt
merene hodnoty = zap.txt
... nebo po stisknutí 2 můžeš zadat vstupní hodnoty (čísla bodů, souřadnice y a x, úhly a délky) ručně
Nacist data ze souboru - 1
Zadat data rucne - 2
Typ vstupu : 2
Pocatecni pripojovaci bod
cislo = 101
y [m] = 726982.256
x [m] = 1026352.842
Orientace na bod
cislo =
Po zadání vstupních dat (ať již ze souboru, nebo ručně) je ještě třeba zvolit, zda se budou počítat apriorní nebo aposteriorní střední chyby a parametry středních elips chyb
Vypocet APRIORNICH strednich chyb a parametru elips - 1
Vypocet APOSTERIORNICH strednich chyb a parametru elips - 2
Typ vypoctu :
Poté je proveden vlastní výpočet (objeví se při něm výzva k zadání apriorní střední chyby jednotkové m0, střední chyby délek md a úhlů mw) jehož výsledkem je níže uvedený protokol a obrázek polygonu se středními elipsami chyb
Příklad
Vstupní soubory:
seznam souřadnic (číslo bodu, souřadnice y, souřadnice x)
Čísla bodů mohou obsahovat čísla i písmena, nesmí tam ale být mezery
Desetinným oddělovačem je tečka
Q 733762.80 1030209.84
1 733590.48 1029780.16
8 733105.50 1028711.45
L 733313.52 1028743.33
zápisník měřených hodnot (číslo bodu, úhel na bodě, délka vpřed)
1 249.4886 397.209
2 244.4452 283.242
3 089.1833 379.817
4 240.3127 267.978
5 150.8930 300.268
6 140.1824 254.409
7 149.7205 164.782
8 201.8150
Uvedené hodnoty jsou z příkladu ve skriptech Geodézie 50 - Příklady a návody na cvičení na straně 86 až 91 (výsledné střední chyby a délky poloos elips chyb se liší, protože v příkladu počítají aposteriorní hodnoty)
Výsledkem výpočtu je takovýto protokol ...
========================================================================
VYROVNANI OBOUSTRANE PRIPOJENEHO A ORIENTOVANEHO POLYGONOVEHO PORADU
========================================================================
apriorni stredni chyba jednotkova m0 [mm/cc] = 15
stredni chyba delek md [mm] = 15
stredni chyba uhlu mw [cc] = 20
Dane a merene hodnoty
---------------------
bod y[m] x[m] uhel[g] delka[m]
-----------------------------------------------------------------
Q 733762.8000 1030209.8400
1 733590.4800 1029780.1600 249.4886 397.209
2 244.4452 283.242
3 89.1833 379.817
4 240.3127 267.978
5 150.8930 300.268
6 140.1824 254.409
7 149.7205 164.782
8 733105.5000 1028711.4500 201.8150
L 733313.5200 1028743.3300
Uzavery (matice podminek):
--------------------------
Ux = -48.59 mm
Uy = +16.03 mm
Uw = +28.48 cc
Kontroly vypoctu
----------------
vTPv = 481.159304238
-UTK = 481.159304238
Aposteriorni stredni chyba jednotkova m0
----------------------------------------
m0 = 12.66 [mm/cc]
Dosazeni do pretvorenych podminkovych rovnic:
---------------------------------------------
Ux = +0.000000000000007
Uy = -0.000000000000011
Uw = +0.000000000000021
Mezni uzavery na hladine vyznamnosti 95%:
-----------------------------------------
uzaver mezni uzaver vysledek
------------------------------------------------------
Ux = -48.59 mm Mx = ±90.52 mm splneno
Uy = +16.03 mm My = ±164.67 mm splneno
Uw = +28.48 cc Mw = ±141.42 cc splneno
Vyrovnane merene hodnoty a jejich apriorni stredni chyby
--------------------------------------------------------
bod uhel[g] m[cc] delka[m] m[mm]
---------------------------------------------------
1 249.48979 14.94 397.2085 13.89
2 244.44554 17.79 283.2465 14.28
3 89.18315 17.39 379.8101 13.45
4 240.31219 18.51 267.9745 13.56
5 150.89198 17.96 300.2607 13.52
6 140.18118 17.52 254.4028 14.22
7 149.71959 17.31 164.7808 14.08
8 201.81444 17.18
Vyrovnane souradnice urcovanych bodu, jejich stredni chyby a parametry elips chyb
---------------------------------------------------------------------------------
bod y[m] x[m] my[mm] mx[mm] mp[mm] mxy[mm] omega[g] a[mm] b[mm]
----------------------------------------------------------------------------------------------------------
2 733226.5105 1029621.0979 13.26 10.20 16.73 11.83 73.6064 13.89 9.32
3 732954.7807 1029701.0449 18.22 14.66 23.39 16.54 77.4553 18.72 14.02
4 732910.7346 1029323.7974 16.66 16.26 23.28 16.46 64.9828 16.90 16.01
5 732728.1801 1029127.6236 16.64 17.75 24.33 17.20 194.4874 17.76 16.63
6 732734.7656 1028827.4351 16.81 13.39 21.49 15.20 128.8071 17.74 12.13
7 732943.3984 1028681.8566 14.00 4.69 14.77 10.44 92.0260 14.10 4.39
... a obrázek s náčrtem polygonového pořadu se středními elipsami chyb

Zdrojový kód programu
function plg(ss,zap)
%
%+------------------------------------------------------------------------+
%| VYROVNANI OBOUSTRANE PRIPOJENEHO A ORIENTOVANEHO POLYGONOVEHO PORADU |
%+------------------------------------------------------------------------+
%
%funkci spust ve tvaru PLG('SEZNAM SOURADNIC' , 'ZAPISNIK MERENYCH HODNOT')
%
%jmena souboru zadat vcetne apostrofu '
%
%Pokud funkci spustis bez parametru, zobrazi se nabidka s moznostmi
%zadat vstupni soubory (seznam souradnic a zapisnik merenych hodnot)
%nebo zadat vstupni hodnoty (cisla bodu, souradnice y a x a merene
%uhly a delky) rucne
%
%Vstupni soubory musi byt ve stejnem adresari jako tato funkce (soubor plg.m),
%jinak je treba zadat uplnou cestu, napriklad PLG('G:/SS.TXT','ZAP.TXT')
%
%Pri prvnim pouziti je asi nejlepsi zadat hodnoty rucne a pote se podivat
%na struktutu souboru 'ss.txt' a 'zap.txt'
%
%
%Napsano a odladeno v Matlabu 6 (R12), nezarucuji tudiz, ze to ve starsich
%verzich bude chodit
%
%
%rsc@email.cz
if (nargin == 2)
%nacteni vstupnich hodnot ze zadanych souboru (pokud existuji)
%------------------------------------------
elseif (nargin ~= 2)
clc;
disp('Nacist data ze souboru - 1');
disp(' Zadat data rucne - 2');
c = input('\n Typ vstupu : ');
if c==1
disp(' ');
disp(' ');
disp('Nacteni dat ze souboru');
disp('(jmeno souboru zadavej ve tvaru jmeno.pripona (bez''))');
ss = input(' seznam souradnic = ','s');
zap = input(' merene hodnoty = ','s');
%------------------------------------------
%rucni zadani vstupnich hodnot
elseif c==2
disp(' ');
disp(' ');
%---------------------------------
%pripojovaci body
disp('Pocatecni pripojovaci bod');
ss_c(2,:) = input(' cislo = ','s');
d = length(ss_c(2,:));
ss(2,1) = input(' y [m] = ');
ss(2,2) = input(' x [m] = ');
%---------------------------------
disp(' ');
disp('Orientace na bod');
cislo = input(' cislo = ','s');
if(length(cislo) > d)
ss_c(1,length(cislo)) = ' ';
d = length(cislo);
end
if(length(cislo) < d)
cislo(1,d) = ' ';
end
ss_c(1,:) = cislo;
ss(1,1) = input(' y [m] = ');
ss(1,2) = input(' x [m] = ');
%---------------------------------
disp(' ');
disp(' ');
disp('Koncovy pripojovaci bod');
cislo = input(' cislo = ','s');
if(length(cislo) > d)
ss_c(3,length(cislo)) = ' ';
d = length(cislo);
end
if(length(cislo) < d)
cislo(1,d) = ' ';
end
ss_c(3,:) = cislo;
ss(3,1) = input(' y [m] = ');
ss(3,2) = input(' x [m] = ');
%---------------------------------
disp(' ');
disp('Orientace na bod');
cislo = input(' cislo = ','s');
if(length(cislo) > d)
ss_c(4,length(cislo)) = ' ';
d = length(cislo);
end
if(length(cislo) < d)
cislo(1,d) = ' ';
end
ss_c(4,:) = cislo;
ss(4,1) = input(' y [m] = ');
ss(4,2) = input(' x [m] = ');
%zapis do souboru "ss.txt"
fid = fopen('ss.txt','w');
for(i=1:4)
fprintf(fid,'%s \t %f \t %f\n',ss_c(i,:),ss(i,1),ss(i,2));
end
fclose(fid);
disp(' ');
disp('Zadane hodnoty byly zapsany do souboru "ss.txt"');
%merene uhly a delky
disp(' ');
disp(' ');
p = input('Pocet vrcholu polygonoveho poradu = ');
%---------------------------------
disp(' ');
disp('Vrchol 1 (pocatecni pripojovaci bod)');
zap_c(1,:) = ss_c(2,:); %cislo prvniho bodu
d = length(zap_c(1,:));
zap(1,1) = input(' uhel [g] = ');
zap(1,2) = input(' delka vpred [m] = ');
%---------------------------------
for i=2:p-1
fprintf(1,'\nVrchol %i\n',i);
cislo = input(' cislo = ','s');
if(length(cislo) > d)
zap_c(1,length(cislo)) = ' ';
d = length(cislo);
end
if(length(cislo) < d)
cislo(1,d) = ' ';
end
zap_c(i,:) = cislo;
zap(i,1) = input(' uhel [g] = ');
zap(i,2) = input(' delka vpred [m] = ');
end
%---------------------------------
fprintf(1,'\nVrchol %i (koncovy pripojovaci bod)\n',p);
cislo = ss_c(3,:); %cislo posledniho bodu
if(length(cislo) > d)
ss_c(1,length(cislo)) = ' ';
d = length(cislo);
end
if(length(cislo) < d)
cislo(1,d) = ' ';
end
zap_c(p,:) = cislo;
zap(p,1) = input(' uhel [g] = ');
%zapis do souboru "zap.txt"
fid = fopen('zap.txt','w');
for(i=1:p-1)
fprintf(fid,'%s \t %f \t %f\n',zap_c(i,:),zap(i,1),zap(i,2));
end
fprintf(fid,'%s \t %f\n',zap_c(p,:),zap(p,1));
fclose(fid);
disp(' ');
disp('Zadane hodnoty byly zapsany do souboru "zap.txt"');
clear ss;
clear zap;
clear ss_c;
clear zap_c;
ss = [];
zap = [];
ss_c = [];
zap_c = [];
end
end
%------------------------------------------
%vlastni nacteni vstupnich dat
if(isempty(ss) == 1) ss = 'ss.txt';
end
if(isempty(zap) == 1) zap = 'zap.txt';
end
%body
[fid,mes]=fopen(ss,'r');
if (fid < 0)
fprintf(2,'PRI NACTENI SEZNAMU SOURADNIC DOSLO K CHYBE\n',ss,mes);
error('')
end
[p1,p2]=fscanf(fid,'%s');
fseek(fid,0,-1);
poc_b = p2/3;
if (rem(p2,3)~=0)
error('Chyba v seznamu souradnic');
end
for i=1:poc_b
pom = fscanf(fid,'%s',1);
[m,n] = size(pom);
cis(i,1:n) = pom;
y_dane(i,1) = fscanf(fid,'%f',1);
x_dane(i,1) = fscanf(fid,'%f',1);
end
%merene hodnoty
[f,mes]=fopen(zap,'r');
if f<0
fprintf(2,'Soubor se nepodarilo otevrit\n',soubor,mes);
error('')
end
[p1,p2]=fscanf(f,'%s');
fseek(f,0,-1);
poc_mereni = ((p2+1)/3);
if (rem((p2+1),3) ~= 0)
error('Chyba v zapisniku merenych dat');
end
for i=1:poc_mereni
pom = fscanf(f,'%s',1);
[m,n] = size(pom);
cis_mer(i,1:n) = pom;
mer_uhly(i,1) = fscanf(f,'%f',1);
if(i ~= poc_mereni)
mer_delky(i,1) = fscanf(f,'%f',1);
end
end
p = poc_mereni; %p = pocet vrcholu poradu
%vektor cisel bodu
for i=1:p
cisla_bodu(i+1,:) = cis_mer(i,:);
end
[dd d] = size(cisla_bodu);
cislo = cis(1,:);
if(length(cislo) > d)
cisla_bodu(1,length(cislo)) = ' ';
end
if(length(cislo) < d)
cislo(1,d) = ' ';
end
cisla_bodu(1,:) = cislo;
cislo = cis(4,:);
if(length(cislo) > d)
cisla_bodu(1,length(cislo)) = ' ';
end
if(length(cislo) < d)
cislo(1,d) = ' ';
end
cisla_bodu(p+2,:) = cislo;
%******************************************************
typ_m0 = 1; %implicitni hodnota
disp(' ');
disp(' Vypocet APRIORNICH strednich chyb a parametru elips - 1');
disp('Vypocet APOSTERIORNICH strednich chyb a parametru elips - 2');
typ_m0 = input('\n Typ vypoctu : ');
%******************************************************
clc;
disp('========================================================================');
disp(' VYROVNANI OBOUSTRANE PRIPOJENEHO A ORIENTOVANEHO POLYGONOVEHO PORADU ');
disp('========================================================================');
disp(' ')
disp(' ')
ro = 200/pi;
roq = 2000/pi;
%apriorni stredni chyby
m0 = input('apriorni stredni chyba jednotkova m0 [mm/cc] = ');
m_d = input(' stredni chyba delek md [mm] = ');
m_w = input(' stredni chyba uhlu mw [cc] = ');
%m0 = 15;
%m_d = 15;
%m_w = 20;
disp(' ')
disp(' ')
disp('Dane a merene hodnoty');
disp('---------------------');
disp('bod y[m] x[m] uhel[g] delka[m]');
disp('-----------------------------------------------------------------');
fprintf(1,'%s \t %.4f \t %.4f\n',cisla_bodu(1,:),y_dane(1,1),x_dane(1,1));
if(mer_uhly(1,1) <= 100)
fprintf(1,'%s \t %.4f \t %.4f \t \t %.4f \t %.3f \n',cisla_bodu(2,:),y_dane(2,1),x_dane(2,1),mer_uhly(1,1),mer_delky(1,1));
end
if(mer_uhly(1,1) >= 100)
fprintf(1,'%s \t %.4f \t %.4f \t \t %.4f \t %.3f \n',cisla_bodu(2,:),y_dane(2,1),x_dane(2,1),mer_uhly(1,1),mer_delky(1,1));
end
for i=2:p-1
if(mer_uhly(i,1) <= 100)
fprintf(1,'%s \t \t \t \t \t \t \t \t %.4f \t %.3f\n',cisla_bodu(i+1,:),mer_uhly(i,1),mer_delky(i,1));
end
if(mer_uhly(i,1) >= 100)
fprintf(1,'%s \t \t \t \t \t \t \t \t %.4f \t %.3f\n',cisla_bodu(i+1,:),mer_uhly(i,1),mer_delky(i,1));
end
end
if(mer_uhly(p,1) <= 100)
fprintf(1,'%s \t %.4f \t %.4f \t \t %.4f\n',cisla_bodu(p+1,:),y_dane(3,1),x_dane(3,1),mer_uhly(p,1));
end
if(mer_uhly(p,1) >= 100)
fprintf(1,'%s \t %.4f \t %.4f \t \t %.4f\n',cisla_bodu(p+1,:),y_dane(3,1),x_dane(3,1),mer_uhly(p,1));
end
fprintf(1,'%s \t %.4f \t %.4f\n',cisla_bodu(p+2,:),y_dane(4,1),x_dane(4,1));
disp(' ');
disp(' ');
disp(' ');
disp('Uzavery (matice podminek):');
disp('--------------------------');
%smernik PQ
sigma_pq = ro*atan2(y_dane(1,1) - y_dane(2,1),x_dane(1,1) - x_dane(2,1));
if(sigma_pq < 0)
sigma_pq = sigma_pq + 400;
end
%smernik KL
sigma_kl = ro*atan2(y_dane(4,1) - y_dane(3,1),x_dane(4,1) - x_dane(3,1));
if(sigma_kl < 0)
sigma_kl = sigma_kl + 400;
end
%priblizne smerniky a souradnicove rozdily
pribl_smer(1,1) = sigma_pq + mer_uhly(1,1);
for i=2:p
pribl_y(i-1,1) = mer_delky(i-1,1) * sin(pribl_smer(i-1,1)/ro);
pribl_x(i-1,1) = mer_delky(i-1,1) * cos(pribl_smer(i-1,1)/ro);
pribl_smer(i,1) = pribl_smer(i-1,1) + mer_uhly(i,1) - 200;
end
%uzaver sour x
U_x = 1000*(x_dane(2,1) + sum(pribl_x) - x_dane(3,1));
U(1,1) = U_x;
fprintf(1,'\tUx = %+.2f mm\n', U_x);
%uzaver sour y
U_y = 1000*(y_dane(2,1) + sum(pribl_y) - y_dane(3,1));
U(2,1) = U_y;
fprintf(1,'\tUy = %+.2f mm\n', U_y);
%uhlovy uzaver
U_w = 10000*(sigma_pq + sum(mer_uhly) - (200*round((sigma_pq + sum(mer_uhly) - sigma_kl)/200)) - sigma_kl);
U(3,1) = U_w;
fprintf(1,'\tUw = %+.2f cc\n\n\n', U_w);
%matice D
%delky
for i=1:p-1
D(1,i) = cos(pribl_smer(i,1)/ro);
D(2,i) = sin(pribl_smer(i,1)/ro);
end
%uhly
dy = sum(pribl_y);
dx = sum(pribl_x);
for i=p:2*p-2
for j=1:i-p+1
D(1,i) = -dy/roq;
D(2,i) = dx/roq;
D(3,i) = 1;
end
dy = dy - pribl_y(j,1);
dx = dx - pribl_x(j,1);
end
D(3,2*p-1) = 1;
%matice vah P
%delky
for i=1:p-1
P(i,i) = (m0^2)/(m_d^2);
end
%uhly
for i=p:2*p-1
P(i,i) = (m0^2)/(m_w^2);
end
%matice normalnich rovnic
N = D*inv(P)*D';
%korelaty
K = -inv(N)*U;
%vektor oprav merenzch hodnot
v = inv(P)*D'*K;
%kontroly
disp(' ');
disp('Kontroly vypoctu');
disp('----------------');
%vTPv
vTPv = v'*P*v;
fprintf(1,'\tvTPv = %.9f \n', vTPv);
%-UTK
UTK = -U'*K;
fprintf(1,'\t-UTK = %.9f \n\n', UTK);
%aposteriorni m0
disp('Aposteriorni stredni chyba jednotkova m0');
disp('----------------------------------------');
apost_m0 = (vTPv/3)^0.5;
fprintf(1,'\t m0 = %.2f [mm/cc]\n\n', apost_m0);
if(typ_m0 == 2) m0 = apost_m0;
end
%zpetne dosazeni do PPR (=0)
disp('Dosazeni do pretvorenych podminkovych rovnic:');
disp('---------------------------------------------');
PPR = D*v+U;
fprintf(1,'\t Ux = %+.15f \n', PPR(1,1));
fprintf(1,'\t Uy = %+.15f \n', PPR(2,1));
fprintf(1,'\t Uw = %+.15f \n\n\n', PPR(3,1));
%vyrovnane delky
for i=1:p-1
vyr_delky(i,1) = mer_delky(i,1) + v(i,1)/1000;
end
%vyrovnane uhly
for i=1:p
vyr_uhly(i,1) = mer_uhly(i,1) + v(i+p-1,1)/10000;
end
%mezni uzavery
disp('Mezni uzavery na hladine vyznamnosti 95%:');
disp('-----------------------------------------');
Mu = 2.5*sqrt(m0^2*N);
disp(' uzaver mezni uzaver vysledek');
disp('------------------------------------------------------');
if(U_x <= Mu(1,1))
fprintf(1,'\t Ux = %+.2f mm \t Mx = ±%.2f mm \t splneno\n',U_x,Mu(1,1));
end
if(U_x > Mu(1,1))
fprintf(1,'\t Ux = %+.2f mm \t Mx = ±%.2f mm \t NESPLNENO\n',U_x,Mu(1,1));
end
if(U_y <= Mu(2,2))
fprintf(1,'\t Uy = %+.2f mm \t My = ±%.2f mm \t splneno\n',U_y,Mu(2,2));
end
if(U_y > Mu(2,2))
fprintf(1,'\t Uy = %+.2f mm \t My = ±%.2f mm \t NESPLNENO\n',U_y,Mu(2,2));
end
if(U_w <= Mu(3,3))
fprintf(1,'\t Uw = %+.2f cc \t Mw = ±%.2f cc \t splneno\n\n',U_w,Mu(3,3));
end
if(U_w > Mu(3,3))
fprintf(1,'\t Uw = %+.2f cc \t Mw = ±%.2f cc \t NESPLNENO\n\n',U_w,Mu(3,3));
end
%kovariancni matice vyrovnanych mereni
Qt = inv(P)-inv(P)*D'*inv(N)*D*inv(P);
Mt = (m0^2)*Qt;
%vypis
if(typ_m0 == 2)
disp('Vyrovnane merene hodnoty a jejich aposteriorni stredni chyby');
disp('------------------------------------------------------------');
end
if (typ_m0 ~= 2)
disp('Vyrovnane merene hodnoty a jejich apriorni stredni chyby');
disp('--------------------------------------------------------');
end
disp('bod uhel[g] m[cc] delka[m] m[mm]');
disp('---------------------------------------------------');
for i=1:p-1
fprintf(1,'%s \t %.5f \t %.2f \t %.4f \t %.2f\n',cisla_bodu(i+1,:),vyr_uhly(i,1),(Mt(i+p-1,i+p-1)^0.5),vyr_delky(i,1),(Mt(i,i))^0.5);
end
fprintf(1,'%s \t %.5f \t %.2f\n\n\n',cisla_bodu(p+1,:),vyr_uhly(p,1),(Mt(2*p-1,2*p-1)^0.5));
%vyrovnane smerniky a souradnicove rozdily
vyr_smer(1,1) = sigma_pq + vyr_uhly(1,1);
for i=2:p
vyr_smer(i,1) = vyr_smer(i-1,1) + vyr_uhly(i,1) - 200;
vyr_y(i-1,1) = vyr_delky(i-1,1) * sin(vyr_smer(i-1,1)/ro);
vyr_x(i-1,1) = vyr_delky(i-1,1) * cos(vyr_smer(i-1,1)/ro);
end
%vyrovnane souradnice urcovanych bodu
y_vyr(1,1) = y_dane(2,1) + vyr_y(1,1);
x_vyr(1,1) = x_dane(2,1) + vyr_x(1,1);
for i=2:p-2
y_vyr(i,1) = y_vyr(i-1,1) + vyr_y(i,1);
x_vyr(i,1) = x_vyr(i-1,1) + vyr_x(i,1);
end
%matice F
%delky
for i=1:p-2
for j=i:p-2
F(2*j-1,i) = cos(vyr_smer(i,1)/ro);
F(2*j,i) = sin(vyr_smer(i,1)/ro);
end
end
%uhly
for i=1:p-2
dy = 0;
dx = 0;
for j=i:p-2
dy = dy + vyr_y(j,1);
F(2*j-1,i+p-1) = -dy/roq;
dx = dx + vyr_x(j,1);
F(2*j,i+p-1) = dx/roq;
end
end
F(3,2*p-1) = 0;
%kovariancni matice souradnic urcovanych bodu
Mx = F*Mt*F';
%elipsy chyb
d=diag(Mx);
mx=d(1:2:end);
my=d(2:2:end);
d2=diag(Mx,1);
mxy=d2(1:2:end);
w=0.5*atan2((2*(mxy)),(mx-my));
a=sqrt(mx.*cos(w).^2+mxy.*sin(2.*w)+my.*sin(w).^2);
b=sqrt(mx.*sin(w).^2-mxy.*sin(2.*w)+my.*cos(w).^2);
for i=1:p-2
w_g(i,1) = ro*w(i,1);
if(w(i,1) < 0)
w_g(i,1) = 200+ro*w(i,1);
end
end
%vypis
if(typ_m0 == 2)
disp('Vyrovnane souradnice urcovanych bodu, jejich aposteriorni stredni chyby a parametry elips chyb');
disp('----------------------------------------------------------------------------------------------');
end
if (typ_m0 ~= 2)
disp('Vyrovnane souradnice urcovanych bodu, jejich apriorni stredni chyby a parametry elips chyb');
disp('------------------------------------------------------------------------------------------');
end
disp('bod y[m] x[m] my[mm] mx[mm] mp[mm] mxy[mm] omega[g] a[mm] b[mm]');
disp('----------------------------------------------------------------------------------------------------------');
for i=1:p-2
fprintf(1,'%s \t %.4f \t %.4f \t %.2f \t %.2f \t %.2f \t %.2f \t %.4f \t %.2f \t %.2f \n',cisla_bodu(i+2,:),y_vyr(i,1),x_vyr(i,1),Mx(2*i,2*i)^0.5,Mx(2*i-1,2*i-1)^0.5,(Mx(2*i,2*i) + Mx(2*i-1,2*i-1))^0.5,((Mx(2*i,2*i) + Mx(2*i-1,2*i-1))^0.5)/(2^0.5),w_g(i,1),a(i,1),b(i,1));
end
disp(' ');
disp(' ');
disp(' ');
disp(' ');
disp(' ');
%graf
%body
x_P(1,1) = x_dane(1,1);
x_P(2,1) = x_dane(2,1);
y_P(1,1) = y_dane(1,1);
y_P(2,1) = y_dane(2,1);
y_T(1,1) = y_dane(2,1);
x_T(1,1) = x_dane(2,1);
for i=1:p-2
x_T(i+1,1) = x_vyr(i);
y_T(i+1,1) = y_vyr(i);
end
y_T(p,1) = y_dane(3,1);
x_T(p,1) = x_dane(3,1);
x_K(1,1) = x_dane(3,1);
x_K(2,1) = x_dane(4,1);
y_K(1,1) = y_dane(3,1);
y_K(2,1) = y_dane(4,1);
hold on;
axis equal;
%polygon
plot(-y_P,-x_P,'b:');
plot(-y_T,-x_T,'k-.');
plot(-y_K,-x_K,'b:');
%popis bodu
x = [];
y = [];
x(1,1) = x_dane(1,1);
x(2,1) = x_dane(2,1);
y(1,1) = y_dane(1,1);
y(2,1) = y_dane(2,1);
for i=1:p-2
x(i+2,1) = x_vyr(i);
y(i+2,1) = y_vyr(i);
end
x(p+1,1) = x_dane(3,1);
x(p+2,1) = x_dane(4,1);
y(p+1,1) = y_dane(3,1);
y(p+2,1) = y_dane(4,1);
for i=1:p+2
text(-y(i,1),-x(i,1),cisla_bodu(i,:),'color','b','spanWeight','bold','spanAngle','italic')
end
x = [];
y = [];
%elipsy
%zvetseni eli
M = ((sum(mer_delky)/(p-1))/2)/((max(a)+max(b))/2);
%poloosy v meritku
a_M = M*a;
b_M = M*b;
w_M = pi-w; %osy GEO -> MAT
%vykresleni
for i=1:p-2
t=0:0.1:2*pi+0.1;
x1 = a_M(i,1) * cos(t);
y1 = b_M(i,1) * sin(t);
x2 = x1 * cos(w_M(i,1)) + y1 * sin(w_M(i,1)) + x_vyr(i,1);
y2 = -x1 * sin(w_M(i,1)) + y1 * cos(w_M(i,1)) + y_vyr(i,1);
%vykresleni eli
plot(-y2,-x2,'r');
%osy eli -> a
xA = x_vyr(i,1) + a_M(i,1) * cos(w(i,1));
xB = x_vyr(i,1) - a_M(i,1) * cos(w(i,1));
yA = y_vyr(i,1) + a_M(i,1) * sin(w(i,1));
yB = y_vyr(i,1) - a_M(i,1) * sin(w(i,1));
y=[yA,yB];
x=[xA,xB];
plot(-y,-x,'r:');
%osy eli -> b
xA = x_vyr(i,1) + b_M(i,1) * cos(w(i,1) + pi/2);
xB = x_vyr(i,1) - b_M(i,1) * cos(w(i,1) + pi/2);
yA = y_vyr(i,1) + b_M(i,1) * sin(w(i,1) + pi/2);
yB = y_vyr(i,1) - b_M(i,1) * sin(w(i,1) + pi/2);
y=[yA,yB];
x=[xA,xB];
plot(-y,-x,'r:');
end
title('Polygonový porad a elipsy chyb','spanWeight','bold','spanSize',12,'spanAngle','italic','VerticalAlignment','bottom');
axis off;
hold off;