function amplappr(LTIsys)
% function amplappr(LTIsys)
%
% Zeichnet die Geradenapproximation von
% Übertragungsfunktionen.
%
% J.Schröder, J.Lunze
% für 3. Auflage: 14.2.2001
% für 4. Auflage: 15.9.2003
% für 5. Auflage: 18.4.2005
% für 6. Auflage: 14.12.2006
% für 9. Auflage: 13.2.2012
% für 12. Auflage: 03.07.2019

LTItf=tf(LTIsys);
num=LTItf.num{:};
den=LTItf.den{:};
zpk_mod=zpk(LTItf);
k_stat=zpk_mod.k;

zeros=roots(num);
poles=roots(den);

%   fh=figure;
%   set(fh,'Name','Bode Diagramm, Geradenapproximation');
%   set(fh,'NumberTitle','off');
%   clf

if ~isempty(zeros)
   zt=zeros+max([1 max(abs(zeros))])*(abs(zeros)==0);
else
   zt=[];
end
if ~isempty(poles)
   pt=poles+max([1 max(abs(poles))])*(abs(poles)==0);
else
   pt=[];
end

wmin=min([abs(zt); abs(pt)])/10;
wmax=max([abs(zeros); abs(poles)])*10;

% Kuerzen
z=zeros;
p=poles;
LTItf=tf(zpk(zeros,poles,k_stat));
num=LTItf.num{:};
den=LTItf.den{:};

for k=1:length(zeros)
  if sum(p==zeros(k))>0
    indp=find(p==zeros(k));
    if ~isempty(z)
      indz=find(z==zeros(k));
    end
    if ~isempty(indz)
      z=cutout(z,indz(1:min(length(indz),length(indp))));
    end
    if ~isempty(indp)
      p=cutout(p,indp(1:min(length(indz),length(indp))));
    end
  end
end
zeros=z;
poles=p;

wknick=[sort(abs([zeros; poles])); wmax];

if isempty(zeros)
  d_part=0;
else
  d_part=sum(abs(zeros)==0);
end
if isempty(poles)
  i_part=0;
else
  i_part=sum(abs(poles)==0);
end

steigung=20*(d_part-i_part);

w=[wmin];
ampl=[];
phas=[];
knickcount=1;
if steigung~=0
   w_0dB=(1/k_stat)^(1/(d_part-i_part));
   % Geradengleichung: y = steigung*log(x) + a
   %               --> log(x) = -a/steigung
   %               --> -a/steigung = log(w_0dB)
   %               --> a = -log(w_0dB)*steigung
   knickcount=knickcount+d_part+i_part;
   abschnitt=-log10(w_0dB)*steigung;
   ampl=[steigung*log10(wmin)+abschnitt,...
steigung*log10(wknick(knickcount))+abschnitt];
   phas=[(d_part-i_part)*90 (d_part-i_part)*90];
   w=[w wknick(knickcount)];
else
   ampl=[20*log10(k_stat) 20*log10(k_stat)];
   phas=[0 0];
   w=[w wknick(knickcount)];
end

while knickcount<length(wknick)
   wk=wknick(knickcount);
   if isempty(zeros)
     nz=0;
   else
     nz=sum(wk==abs(zeros));
   end
   np=sum(wk==abs(poles));
   if nz==0
      steigung=steigung-20*np;
      w=[w wknick(knickcount+np)];
      ampl=[ampl,...
ampl(length(ampl))+steigung*(log10(wknick(knickcount+np))-log10(wk))];
      phas=[phas phas(length(phas))-90*np phas(length(phas))-90*np];
      knickcount=knickcount+np;
   else
      steigung=steigung+20*nz;
      w=[w wknick(knickcount+nz)];
      ampl=[ampl,...
ampl(length(ampl))+steigung*(log10(wknick(knickcount+nz))-log10(wk))];
      phas=[phas phas(length(phas))+90*np phas(length(phas))+90*nz];
      knickcount=knickcount+nz;
   end
end

%subplot(2,1,1);
%cla
semilogx(w,ampl,'b');
a=axis;
a(3)=a(3)-sign(a(3))*0.2*a(3);
a(4)=a(4)+sign(a(4))*0.2*a(4);
axis(a);
latextitle('Geradenapproximation');
grid on;
latexylabel('$$|G|$$ in dB')
latexxlabel('$$\omega$$ in rad/s$$')


%%%%%%%%%%%%%%% SUBFUNS %%%%%%%%%%%%%%%
function vres=cutout(vek,ind)

if length(vek)==0
  vres=[];
  return;
end
v1=[];
if ind>1
  v1=vek(1:ind1-1);
end
v2=[];
if ind<length(vek)
  v2=vek(ind+1:length(vek));
end
vres=[v1; v2];
