echo off; clear; format compact;
%TCV I, 2006/2007, uloha c. 4 
%Zadani c.102 
 
%Priklad 1: 
 
%Odkrokovane[m]     Merene[m] 
x = [ 
     12               11.73 
     16               15.73 
     20               19.78 
] 
%Priklad 2: 
 
alfa=    146.2191; %g
omega=   143.1415; %g
d_AB=      73.773; %m
d_BC=     153.058; %m
m_d=        0.005; %m
m_omega=   0.0011; %g
m_alfa=      0.11; %g

%RESENI
%==========================================================================
%==========================================================================
'Priklad 1:'

%Skutecne chyby krokovanych delek
epsilon_l = x(:,2) - x(:,1)

%Pravy uhel urceny z krokovanych delek
R_krokovano = acos ( (x(1,1)^2+x(2,1)^2-x(3,1)^2) / (2*x(1,1)*x(2,1)) ) * 200 / pi

%Pravy uhel urceny z merenych delek
R_mereno =    acos ( (x(1,2)^2+x(2,2)^2-x(3,2)^2) / (2*x(1,2)*x(2,2)) ) * 200 / pi

%Skutecna chyba vytyceneho praveho uhlu
epsilon_y = R_mereno - R_krokovano

%Vektor parcialnich derivaci funkce (zde R_...) podle jendotlivych stran
syms X
A(1,1) = subs( diff(acos ( (X^2+x(2,1)^2-x(3,1)^2) / (2*X*x(2,1)) )),   X, x(1,1));
A(1,2) = subs( diff(acos ( (x(1,1)^2+X^2-x(3,1)^2) / (2*x(1,1)*X) )),   X, x(2,1));
A(1,3) = subs( diff(acos ( (x(1,1)^2+x(2,1)^2-X^2) / (2*x(1,1)*x(2,1)) )),   X, x(3,1));
A

%Priblizny vypocet sutecne chyby vytyceneho praveho uhlu z linearizovane funkce (zde R_...)
epsilon_y_lin = A * epsilon_l * 200 / pi

%==========================================================================
%==========================================================================
pause;
'Priklad 2.1:'

%Zavedení soøadnic bodu A a prevod uhlu na obloukovou míru
yA = 0;
xA = 0;
rad = pi / 200;
alfa =  alfa * rad;
omega = omega * rad;
m_omega = m_omega * rad;
m_alfa = m_alfa * rad;

%Souradnice bodu B a C
yB = d_AB * sin (alfa)
xB = d_AB * cos (alfa)

yC = yB + d_BC * sin ( alfa + omega - pi )
xC = xB + d_BC * cos ( alfa + omega - pi )

%Pro ilustraci
X = [ alfa
      omega
      d_AB
      d_BC ];   
Y = [ yB
      xB
      yC
      xC];

%Kovariancni matice sigma XX
Sigma_XX = [ m_alfa^2 0         0     0
             0        m_omega^2 0     0
             0        0         m_d^2 0
             0        0         0     m_d^2 ];

%Matice parcialnich derivaci A
A = zeros(4);

syms x
A(1,1) = subs( diff(d_AB * sin (x)),   x, alfa);
A(1,3) = subs( diff(x * sin (alfa)),   x, d_AB);

A(2,1) = subs( diff(d_AB * cos (x)),   x, alfa);
A(2,3) = subs( diff(x * cos (alfa)),   x, d_AB);

A(3,1) = subs( diff(d_AB * sin (x) + d_BC * sin ( x + omega - pi )),   x, alfa);
A(3,2) = subs( diff(d_AB * sin (alfa) + d_BC * sin ( alfa + x - pi )),   x, omega);
A(3,3) = subs( diff(x * sin (alfa) + d_BC * sin ( alfa + omega - pi )),   x, d_AB);
A(3,4) = subs( diff(d_AB * sin (alfa) + x * sin ( alfa + omega - pi )),   x, d_BC);

A(4,1) = subs( diff(d_AB * cos (x) + d_BC * cos ( x + omega - pi )),   x, alfa);
A(4,2) = subs( diff(d_AB * cos (alfa) + d_BC * cos ( alfa + x - pi )),   x, omega);
A(4,3) = subs( diff(x * cos (alfa) + d_BC * cos ( alfa + omega - pi )),   x, d_AB);
A(4,4) = subs( diff(d_AB * cos (alfa) + x * cos ( alfa + omega - pi )),   x, d_BC);

%Kovariancni matice souradnic bodu B a C
Sigma_BC = A * Sigma_XX * A'

m_yB = sqrt ( A(1,:) * Sigma_XX * A(1,:)' )
m_xB = sqrt ( A(2,:) * Sigma_XX * A(2,:)' )
m_yC = sqrt ( A(3,:) * Sigma_XX * A(3,:)' )
m_xC = sqrt ( A(4,:) * Sigma_XX * A(4,:)' )

%==========================================================================
%==========================================================================
'Priklad 2.2:'

d_AC = sqrt ( (yC - yA)^2 + (xC - xA)^2 )

%==========================================================================
%==========================================================================
'Priklad 2.3:'

%==========================================================================
disp(' ')
disp('Priklad 2.3.a: (korektne)')

%Pro ilustraci (Na smìrníku alfa délka nazávisí)
X = [ omega
      d_AB
      d_BC ];
  
Sigma_XX = [ m_omega^2 0     0
             0         m_d^2 0
             0         0     m_d^2 ];
  
Y = [ d_AC ];
%  sqrt(  d_AB^2 + d_BC^2 - 2 * d_AB * d_BC * cos(omega)  )

%Vektor derivaci
A_d_AC_a(1,1) = subs( diff( ...
            ...
            sqrt(  d_AB^2 + d_BC^2 - 2 * d_AB * d_BC * cos(x)  ) ...
            ...
            ), x, omega);

A_d_AC_a(1,2) = subs( diff( ...
            ...
            sqrt(  x^2 + d_BC^2 - 2 * x * d_BC * cos(omega)  ) ...
            ...
            ), x, d_AB);
        
A_d_AC_a(1,3) = subs( diff( ...
            ...
            sqrt(  d_AB^2 + x^2 - 2 * d_AB * x * cos(omega)  ) ...
            ...
            ), x, d_BC);
        
%Kovariancni matice (Variance) d_AC
Sigma_d_AC_a = A_d_AC_a * Sigma_XX * A_d_AC_a';
m_d_AC_a = sqrt (Sigma_d_AC_a)


%==========================================================================
%Výpoèet m_sigma_BC

%Smìrník strany BC
sigma_BC = (alfa + omega - pi) / rad

%Pro ilustraci
X = [ alfa
      omega ];

Sigma_alfa_omega = [ m_alfa^2 0
                     0        m_omega^2 ];

Y = [ sigma_BC ];

%Vektor derivací
A_sigma_BC(1,1) = subs( diff(x + omega - pi),   x, alfa);
A_sigma_BC(1,2) = subs( diff(alfa + x - pi),   x, omega);

%Kovariancni matice (Variance) sigma_BC
Sigma_sigma_BC = A_sigma_BC * Sigma_alfa_omega * A_sigma_BC';
m_sigma_BC = sqrt (Sigma_sigma_BC) / rad %v gonech

%==========================================================================
disp(' ')
disp('Priklad 2.3.b: (nekorektne)')
% d_AC

%Pro ilustraci
X = [ yC
      xC ];

Sigma_yC_xC_d_AC_b = [ m_yC^2 0
                       0      m_xC^2 ];

Y = [ d_AC ];

%Vektor derivaci
A_d_AC_b(1,1) = subs( diff(sqrt ( (x - yA)^2 + (xC - xA)^2 )),   x, yC);
A_d_AC_b(1,2) = subs( diff(sqrt ( (yC - yA)^2 + (x - xA)^2 )),   x, xC);

%Kovariancni matice (Variance) d_AC
Sigma_d_AC_b = A_d_AC_b * Sigma_yC_xC_d_AC_b * A_d_AC_b';
m_d_AC_b = sqrt (Sigma_d_AC_b)

%==========================================================================
% sigma_BC
X = [ m_yB
      m_xB
      m_yC
      m_xC ];
  
Y = [ sigma_BC ];

A_sigma_BC_b(1,1) = subs( diff( ...
            ...
            atan(  (yC - x) / (xC - xB)  ) ...
            ...
            ), x, yB);
A_sigma_BC_b(1,2) = subs( diff( ...
            ...
            atan(  (yC - yB) / (xC - x)  ) ...
            ...
            ), x, xB);
A_sigma_BC_b(1,3) = subs( diff( ...
            ...
            atan(  (x - yB) / (xC - xB)  ) ...
            ...
            ), x, yC);
A_sigma_BC_b(1,4) = subs( diff( ...
            ...
            atan(  (yC - yB) / (x - xB)  ) ...
            ...
            ), x, xC);

Sigma_yC_xC_sigma_BC_b = diag(diag(Sigma_BC));

Sigma_sigma_BC_b = A_sigma_BC_b * Sigma_yC_xC_sigma_BC_b * A_sigma_BC_b';
m_sigma_BC_b = sqrt (Sigma_sigma_BC_b) / rad %v gonech
%==========================================================================

disp(' ')
disp('Priklad 2.3.c: (korektne)')

%Pro ilustraci
X = [ yC
      xC ];

Sigma_yC_xC = [ m_yC^2          Sigma_BC(3,4)
                Sigma_BC(4,3)   m_xC^2          ];

Y = [ d_AC ];

%Kovariancni matice (Variance) d_AC
Sigma_d_AC_c = A_d_AC_b * Sigma_yC_xC * A_d_AC_b';
m_d_AC_c = sqrt (Sigma_d_AC_c)


%==========================================================================
% sigma_BC
X = [ m_yB
      m_xB
      m_yC
      m_xC ];
  
Y = [ sigma_BC ];

Sigma_sigma_BC_b = A_sigma_BC_b * Sigma_BC * A_sigma_BC_b';
m_sigma_BC_b = sqrt (Sigma_sigma_BC_b) / rad %v gonech

%==========================================================================
%==========================================================================
'Priklad 3.1:'

% d_AC = sqrt ( (yC - yA)^2 + (xC - xA)^2 )

Sigma_d_AC_3 = 0.02;
A_d_AC_3 = A_d_AC_a;

Sigma_XX_3 = repmat(Sigma_d_AC_3^2/3,1,3) ./ A_d_AC_3.^2;
m_omega_3 = sqrt(Sigma_XX_3(1,1))/rad
m_d_AB_3 = sqrt(Sigma_XX_3(1,2))
m_d_BC_3 = sqrt(Sigma_XX_3(1,3))

%==========================================================================
'Priklad 3.2:'

P = 0.5 * d_AB * d_BC * sin (omega);

m_P = 0.5;

%Vektor derivací
A_P(1,1) = subs( diff(0.5 * d_AB * d_BC * sin (x)),   x, omega);
A_P(1,2) = subs( diff(0.5 * x * d_BC * sin (omega)),   x, d_AB);
A_P(1,3) = subs( diff(0.5 * d_AB * x * sin (omega)),   x, d_BC);

Sigma_XX_P = repmat(m_P^2/3,1,3) ./ A_P.^2;
m_omega_P = sqrt(Sigma_XX_P(1,1))/rad
m_d_AB_P = sqrt(Sigma_XX_P(1,2))
m_d_BC_P = sqrt(Sigma_XX_P(1,3))

%Kontrola
Sigma_XX = [ Sigma_XX_P(1,1) 0 0
             0 Sigma_XX_P(1,2) 0
             0 0 Sigma_XX_P(1,3) ];
%Kovariancni matice (Variance) d_AC
Sigma_P = A_P * Sigma_XX * A_P';
m_P_kontrola = sqrt (Sigma_P)


