function [am, wa,mag,pha,w]=asy5(cit,jmen,w)
% funkce pro kresleni asymptotickeho prubehu, 
% !!! upraveno z verze pro Matlab 4 do verze 5
% !!! kreslici funkce nejsou dokonceny, pouze pro ziskani dat na kresleni prubehu
% vstupem je citatel a jmenovatel soustavy, popripade i pole frekvencni
% ve kterem chceme pracovat (vsechny poly a nuly musi byt uvnitr hodnot pole
% vystupem je vektor amplitud (v dB) a frekvenci v bodech zlomu a v krajnich
% bodech pole
% vicenasobne koreny jsou uvedeny vicenasobne
% dale vraci prubehy amplitudy, faze a frekvence (jak je generuje bode)
% [zlomy_ampl zlomy_f amplituda faze f] = 
%                asy([100 20 1],conv([1 0.2 1],[0.2 1 0]),w)
% Jestlize nejsou pozadovany vystupni hodnoty, tiskne bode + asymptoty
%
% Pozn.: krajni body pole frekvenci je vhodne volit cca o 1 dekadu od krajnich zlomu
%        V pripade, ze v grafu chceme vice prubehu, je potrebne provest
%	 subplot(2,1,1);hold on;subplot(2,1,2);hold on;
%        zobrazeni prubehu bode:   semilogx(w,20*log10(mag))
%                    a asumptot    semilogx(wa,am) do spolecneho grafu

% nejdrive zjistim zesileni
% a casove konstanty citatele a jmenovatele
[zz pp k] = tf2zp(cit,jmen);
% normalizovany tvar prenosove funkce - nejvyssi mocnina 
% citatele ma koeficient jedna, zesileni je v citatli
[num,den] = zp2tf(zz,pp,k);
% zesileni je takove cislo, ze posledni nenulove koeficienty 
% citatele a jmenovatele jsou rovny jedne
for i=length(num):-1:1
  if num(i)~=0 
    kn = num(i);
    break;
  end
end
for i=length(den):-1:1
  if den(i)~=0 
    kc = den(i);
    break;
  end
end
% tohle je zesileni soustavy
kv = kn / kc;
% ted ty poly a nuly
% nejdrive poly - ty jsou v pp
% poly v nule
pn = 0;
for i=1:length(pp)
  if pp(i)==0
    pn = pn + 1;
  else
    break
  end
end
pn;
%uprava pro fazove neminimalni systemy
zz=-abs(zz);
pp=-abs(pp);
%nuly v nule 
nn = 0;
for i=1:length(zz)
  if zz(i)==0
    nn = nn + 1;
  else
    break
  end
end
nn; 
% musim prohledat ostatni poly, zda nejsou komplexni
% pokud jsou, musim vypitvat casovou konstantu
pp;
for i=pn+1:length(pp)
  if imag(pp(i)) ~= 0
% ma imaginarni cast ale jsou komplexne zdruzene 
% -> bereme jen jednou    
    if imag(pp(i))>0
% ksi / T
      ksit = real(pp(i));
      T = -1/sqrt(imag(pp(i))*imag(pp(i))+ksit*ksit);
% dvojnasobny koren v 1/T, predpoklad -> komplexni poly jdou zasebou
% a zacinaji kladnou imaginarni casti
      pp(i) = 1/T;
      pp(i+1) = 1/T ;
    end
  end
end
pp;
% ted komplexni ve jmenovateli
zz;
for i=nn+1:length(zz)
  if imag(zz(i)) ~= 0
% ma imaginarni cast ale jsou komplexne zdruzene 
% -> bereme jen jednou    
    if imag(zz(i))>0
% ksi / T
      ksit = real(zz(i));
      T = -1/sqrt(imag(zz(i))*imag(zz(i))+ksit*ksit);
% dvojnasobny koren v 1/T, predpoklad -> komplexni poly jdou zasebou
% a zacinaji kladnou imaginarni casti
      zz(i) = 1/T;
      zz(i+1) = 1/T ;
    end
  end
end
zz;
% poly a nuly je potreba usporadat vzestupne, kvyli komplexnim korenum
pp = -pp;
zz = -zz;
pp=sort(pp);
zz=sort(zz);
pp = -pp;
zz = -zz;
% pak ziskam jak to kresli bode
if nargin < 3
  [mag,pha,w]=bode(tf(cit,jmen));
else
  [mag,pha,w]=bode(tf(cit,jmen),w);
end
if nargout == 0
%  figure
%  hold off
%  subplot(2,1,2);
%  semilogx(w,pha);
%  grid on
%  subplot(2,1,1);
%  semilogx(w,20*log10(mag))
%  grid on
%  hold on
end
% a vykreslim do bode to co mam
% celkovy pocet zlomu
pocet = length(zz) + length(pp);
% citace aktualniho polu a nuly
cp = pn+1;
cn = nn+1;
% pocitam pocatek zesileni + sklon
amplituda = 20*log10(kv);
% plot([0.01 1],[amplituda amplituda])
sklon = pn-nn;
amplituda = amplituda + log10(1/w(1))*20*sklon;
% vektory pro ulozeni frekvenci zlomu a hodnot asymptot
am=1:pocet+2-(pn+nn);
wa=1:pocet+2-(pn+nn);
citac = 1;
am(citac)=amplituda;
wa(citac)=w(1);
citac = citac+1;
% plot([0.01 1],[amplituda amplituda])
zz;
cn;
%for i=cp:pocet
for i=1:pocet-(pn+nn)
if cn > length(zz)
% uz jen poly
  f=1;
  sklon;
  if wa(citac-1) ~= -pp(cp)
    wa(citac-1);
    -pp(cp)  ;
%    sklon*log10(-pp(cp)/wa(citac-1))*20;
    amplituda;
    amplituda = amplituda - sklon*log10(-pp(cp)/wa(citac-1))*20;
  end
  sklon = sklon + 1;
  am(citac)=amplituda;
  wa(citac)=-pp(cp);
  citac = citac+1;
  cp = cp+1;
elseif cp > length(pp)
% uz jen nuly
   f=2;
   if wa(citac-1) ~= -zz(cn)
    amplituda = amplituda - sklon*log10(-zz(cn)/wa(citac-1))*20;
  end
  sklon = sklon - 1;
  am(citac)=amplituda;
  wa(citac)=-zz(cn);
  citac = citac+1;
  cn = cn + 1; 
elseif zz(cn) > pp(cp)
% na rade je nula
  f=3;
  sklon;
  if wa(citac-1) ~= -zz(cn)
    amplituda = amplituda - sklon*log10(-zz(cn)/wa(citac-1))*20;
  end
  sklon = sklon - 1;
  am(citac)=amplituda;
  wa(citac)=-zz(cn);
  citac = citac+1;
  cn = cn + 1; 
else
% na rade je pol
  f=4;
  if wa(citac-1) ~= -pp(cp)
    amplituda = amplituda - sklon*log10(-pp(cp)/wa(citac-1))*20;
  end
  sklon = sklon + 1;
  am(citac)=amplituda;
  wa(citac)=-pp(cp);
  citac = citac+1;
  cp = cp+1;
end
end
if am(citac-1) ~= w(length(w))
 amplituda = amplituda - sklon*log10(w(length(w))/wa(citac-1))*20;
end
citac;
am(citac)=amplituda;
wa(citac)=w(length(w));
sklon;
wa;
am;
if nargout == 0
  plot (wa,am,'r');
  hold off;
end
mag = mag(:);

