clear all
close all
clc

%Frekvenční osa
freqs = [10 50 100 500 1000 5000 10000 50000 100000];
freqsC = length(freqs);

%Bočníky
Inoms = [0.03 0.1 0.3 1 5 10 20];
Rnoms = [20 6 2 0.6 0.12 0.06 0.03];
phis = [-1.21E-05 -6.60E-06 -1.37E-06 9.63E-07 9.36E-06 1.43E-05 2.83E-05];
fNom = 1000;                        %Frekvence, na které bylo měřeno phi
shC = length(Inoms);

%Kabely
mK = 0.5;                           % délka kabelu
Rk_Core = 40e-3;                    % odpor jádra
Rk_shield = 36e-3;                  % odpor stínění
Rk = mK*(Rk_Core + Rk_shield);      % odpor celkový (na danou délku)
Lk = mK * 1.4e-6;                   % indukčnost kabelu celkový (na danou délku)   
Ck = mK * 100e-12;                  % indukčnost kabelu celkový (na danou délku)   

%Parametry multimetru
Rdmm = 1e7;
Cdmm = 100e-12;

%Parametry zobrazení
initLineWidth = 2;

%Popisky - legendy grafů
legenda = {'Bez vlivu zateze', 'S pridanou kapacitni odporovou zatezi'};
legendCell = eval(['{' sprintf('''Bocnik pro I_{nom} = %5.2f A'' ',Inoms) '}']);

%Alokace paměti
ZmR = zeros(shC, freqsC);
ZmRC = zeros(shC, freqsC); 
ZphiR = zeros(shC, freqsC); 
ZphiRC = zeros(shC, freqsC); 

%Výpočet výstupů modelu
for shI=1:shC %Pro všechny bočníky
    Inom = Inoms(shI);
    Rnom = Rnoms(shI);
    phi = phis(shI);
    
    %Výpočet Cp/Ls z phi
    if (phi<=0) %RC
        Cp = tan(-phi)/2i*fNom*Rnom;
    else %RL
        Ls = tan(phi)/2i*fNom*Rnom;
    end;
    
    for fI = 1:freqsC;  %Pro všechny frekvence
        f = freqs(fI);

        %Výpočet Ak kabelu
        Ak = [1+(Rk + 2i*pi*f*Lk)*2i*pi*f*Ck   Rk + 2i*pi*f*Lk;   ...
              2i*pi*f*Ck                       1                  ];

        %Výpočet Admm multimetru
        AdmmRC = [1 0; 2i*pi*f*Cdmm + 1/Rdmm 1];

        %Výpočet Ash - matice bočníku
        if (phi<=0) %RC
            ASh = [1 0; 1/Rnom + 2i*pi*f*Cp  1];
        else %RL
            ASh = [1 0; Rnom + 2i*pi*f*Ls 1];
        end;
        
        %Celkové matice
        ARC = Ak*ASh*Ak*AdmmRC;
        AR = Ak*ASh*Ak;
        
        %přenosová impedance
        ZcR = 1/(AR(2,1));
        ZcRC = 1/(ARC(2,1));
        % modul impedance
        ZmR(shI, fI) = real(ZcR);
        ZmRC(shI, fI) = real(ZcRC); 
        % fázový posuv
        ZphiR(shI, fI) = atan( imag(ZcR) / real(ZcR) );
        ZphiRC(shI, fI) = atan( imag(ZcRC) / real(ZcRC) );
    end
end

%Převod modulu impedance na mOhm
ZmR_mOhm = ZmR*1e3;
ZmRC_mOhm = ZmRC*1e3;

%Převod fázového posuvu na urad
ZphiR_urad = ZphiR*1e6;
ZphiRC_urad = ZphiRC*1e6;

%Výpočet ABS chyb
ZmAbs = ZmRC - ZmR;
ZphiAbs = ZphiRC - ZphiR;

%Výpočet REL chyb
ZmRel = ZmAbs./ZmR .* 1e6;
ZphiRel = ZphiAbs./ZphiR .* 1e2;

% Úprava rozlišení abs chyby
ZmAbs_nOhm = ZmAbs * 1e6;
ZphiAbs_mikroOhm = ZphiAbs * 1e6;

%Zobrazení
for shI=1:shC
    Inom = Inoms(shI);
    Rnom = Rnoms(shI);
    phi = phis(shI);
    figure('name', ['Bočník pro I_{nom} = ', num2str(Inom)])
    set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
    
    subplot(2, 1, 1);
    semilogx(freqs, ZmR_mOhm(shI, :), '-r', 'LineWidth', initLineWidth);
    hold on;
    semilogx(freqs, ZmRC_mOhm(shI, :), '--b', 'LineWidth', initLineWidth);
    hold off;
    grid on;
    legend(legenda);
    title('Modul impedance');
    xlabel('$f [ Hz ]$', 'interpreter','latex');
    ylabel('$Z_m [ m\Omega ]$', 'interpreter','latex');

    subplot(2, 1, 2);
    semilogx(freqs, ZphiR_urad(shI, :), '-r', 'LineWidth', initLineWidth);
    hold on;
    semilogx(freqs, ZphiRC_urad(shI, :), '--b', 'LineWidth', initLineWidth);
    hold off;
    grid on;
    legend(legenda, 'Location', 'SouthWest');
    title('Fázový posuv');
    xlabel('$f [ Hz ]$', 'interpreter','latex');
    ylabel('$\varphi$ [ $\mu$ rad ]', 'interpreter','latex');
    
    saveas(gcf, ['fig/sh_', num2str(shI)], 'epsc')
end

%Relativní chyby
figRel = figure('name', 'Relativní chyby')
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
subplot(2, 1, 1);
semilogx(freqs, ZmRel, 'LineWidth', initLineWidth);
grid on;
legend(legendCell, 'Location', 'NorthWest');
title('Relativní chyba amplitudy');
xlabel('$f [ Hz ]$', 'interpreter','latex');
ylabel('$\delta Z_m [ ppm ]$', 'interpreter','latex');

subplot(2, 1, 2);
semilogx(freqs, ZphiRel, 'LineWidth', initLineWidth);
grid on;
legend(legendCell, 'Location', 'NorthWest');
title('Relativní chyba fázového posunu');
xlabel('$f [ Hz ]$', 'interpreter','latex');
ylabel('$\delta Z_m [ \% ]$', 'interpreter','latex');

saveas(figRel, 'fig/sh_rel', 'epsc')

%Absolutní chyby
figAbs = figure('name', 'Absolutní chyby');
set(figAbs, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
subplot(2, 1, 1);
semilogx(freqs, ZmAbs_nOhm, 'LineWidth', initLineWidth);
grid on;
legend(legendCell, 'Location', 'SouthEast');
title('Absolutní chyba amplitudy');
xlabel('$f [ Hz ]$', 'interpreter','latex');
ylabel('$\Delta Z_m [ \mu \Omega ]$', 'interpreter','latex');

subplot(2, 1, 2);
semilogx(freqs, ZphiAbs_mikroOhm, 'LineWidth', initLineWidth);
grid on;
legend(legendCell, 'Location', 'SouthWest');
title('Absolutní chyba fázového posunu');
xlabel('$f [ Hz ]$', 'interpreter','latex');
ylabel('$\Delta\varphi [ \mu rad ]$', 'interpreter','latex');

saveas(figAbs, 'fig/sh_abs', 'epsc')