clc
clear all
close all


% parametry kabelu pro CYKY 3x1.5, 3x2.5 a 3x4 
pi = 3.141592;
mio = 4*pi*10^-7; 
miR = 0.99999;
sigma = 58e6; 
a = 0.892e-3; %3x2.5
d = 3.584e-3; 
Er = 4;
Eo = 8.854e-12;

disp('Generátor pøenosové funkce');
disp('Výbìr topologie:');
disp('1 - Malá (5-6 odboèek)');
disp('2 - Støední (7-8 odboèek)');
disp('3 - Velká (9-10 odboèek)');

while true;
    scenar = input('Zvolte topologii (1-3): ');
        if scenar>0 && scenar <4 
            break
        end
end

%---malá topologie (5 - 6 odboèky)-----------------------------------------------
if      scenar==1   
        nah_cislo=rand(1)+5;    
        pocet_odbocek=round(nah_cislo);  
        delka = round((rand(1,pocet_odbocek+1)*2+4)*10)/10;    %delky tras 4-6 m
        dbr = round((rand(1,pocet_odbocek)*2+4)*10)/10;
                  
%---støední topologie (7 - 8 odboèek)-----------------------------------------------               
elseif  scenar==2   
        nah_cislo=(rand(1)+7);    
        pocet_odbocek=round(nah_cislo);  
        delka = round((rand(1,pocet_odbocek+1)*5+5)*10)/10;    %delky tras 5-10 m
        dbr = round((rand(1,pocet_odbocek)*5+5)*10)/10;
        
%---velká topologie (9 - 10 odboèek)----------------------------------------------         
else  
        nah_cislo=(rand(1)+9);    
        pocet_odbocek=round(nah_cislo);
        delka = round((rand(1,pocet_odbocek+1)*10+10)*10)/10;    %delky tras 10-20 m
        dbr = round((rand(1,pocet_odbocek)*10+10)*10)/10;
end    

%Pøiøazení impedance odboèkám

imp= [5 50 150 1000 100000000];
Zbr= zeros(1,pocet_odbocek);

for im=1:1:pocet_odbocek;
    nah=ceil(rand(1)*5);
    Zbr(im)= imp(nah);
    
end    

%impedance zdroje a zátìže
Zs = 50;
ZL = 50;

disp('           ');
disp('Poèet odboèek:');        
disp(pocet_odbocek);
disp('Délky odboèek [m]:');
disp(dbr);
disp('Impedance odboèek Zbr[Ohm]:');
disp(Zbr);
disp('Délky tras [m]:');
disp(delka);

%zvolení frekvenèního rozsahu
minf = 0; % minimálni frekvence
krok = 5000; % krok
maxf = 50e6; % maximálni frekvence
poc = maxf/krok-1;
vysF = zeros(1,poc); % vektor použitých frekvencí
vysH = zeros(1,poc); % vektor výsledné pøenosové funkce

pom = 0; % pomocna promenna

for f= minf:krok:maxf;
    tandelta = (-5.7257*(10^-10*f))+0.06;
    pom = pom + 1;
    omega = 2*pi*f;
    
    % Primarni parametry kabelu 3x2.5
    R = (sqrt(miR*mio.*f./(pi*sigma*a^2)))*((d/(2*a))/sqrt(((d/(2*a))^2)-1));
    L = ((miR*mio)/pi)*acosh(d/(2*a));
    C = (pi*Er*Eo)/(acosh(d/(2*a)));
    G = (2*pi*f*C*tandelta);
    
    % Sekundární parametry
    Zc = sqrt((R+(1i*omega*L))/(G+(1i*omega*C)));
    gama = sqrt((R+(1i*omega*L))*(G+(1i*omega*C)));
   
    % výpoèet ekvivalentní impedace
    Zeq=zeros(1,pocet_odbocek);
    for z=1:1:pocet_odbocek;
        Zeq(z)= Zc*((Zbr(z)+(Zc*tanh(gama*dbr(z))))/(Zc+(Zbr(z)*tanh(gama*dbr(z)))));
    end    

    % výpoèet kaskádních matic
    B = zeros(2,2,pocet_odbocek*2+2);    
    B(:,:,1) = [1 Zs;0 1];    
    i = 1;
    for p=1:1:pocet_odbocek;
        i=i+1;
        B(:,:,i) = [cosh(gama*delka(p)) Zc*sinh(gama*delka(p));(1/Zc)*sinh(gama*delka(p)) cosh(gama*delka(p))];
        i=i+1;
        B(:,:,i) = [1 0;(1/Zeq(p)) 1];
    end   
    i=i+1;
    B(:,:,i) = [cosh(gama*delka(end)) Zc*sinh(gama*delka(end));(1/Zc)*sinh(gama*delka(end)) cosh(gama*delka(end))];
    
    %souèin kaskádních matic
    A=B(:,:,1);
    for j=2:1:i;
        A =A * B(:,:,j);
    end
    
    %výpoèet pøenosové funkce
    H = ZL /((A(1,1)*ZL)+A(1,2)+A(2,1)*ZL+A(2,2)*Zs);    
    vysF(pom) = f; % Vysledny vektor pouzitych frekvenci
    vysH(pom) = H; % Vysledny vektor prenosove fce
    
end

% Vykresleni frekvencni odezvy
plot(vysF/10^6,20*log(abs(vysH))),grid on,title('Frekvenèní odezva'),xlabel('f [MHz]'),ylabel('|H(f)| [dB]');
xlim([0 50]);
legend(sprintf('Model vedení s kaskádními parametry\ngenerovaný pomocí náhodných parametrù'));








