
function [A,vA,w180,wc,GM,PM,Ms,Mt]=minbode_L(num,den,gain,deadtime,w)
%Frequency analysis of SISO system L= gain*deadtime*num/den
%By Truls Larsson, Trondheim 15.09.98 and 21.09.98
%
%Call:  minbode(num,den,gain,deadtime,w) 
%gives:   Bode-plot of L + frequency plot of |S| and |T| (S=1/1+L; T=L/1+L)
%
%or: [A,vA,w180,wc,GM,PM,Ms,Mt]=minbode(num,den,gain,deadtime,w)
%gives:   amplitude (A) og phase (vA) of L as function of frequency,
%       w180 (where va=-180 grader), wc (where A=1), GM, PM, 
%       plus Ms (peak for |S|) and Mt (peak for |T|)
%
%Input:
%num      : numerator polynomial of L
%den      : denominator plynomial of L
%gain     : gain of L
%deadtime : dead time of L
%w        : Frequencies (Note! range must be sufficiently large)
%Example:
%num=[2 1]; den=[25 10 1 0];  gain=0.1; deadtime=3; w= logspace(-2,2,100);
%
if nargin==0
  disp('usage:  minbode(num,den,gain,deadtime,w)')
  disp('or:  [A,vA,w180,wc,GM,PM,Ms,Mt]=minbode(num,den,gain,deadtime,w)')
  return
end

%Frekvensresponsen
teller=polyval(num,w*i);
nevner=polyval(den,w*i);
Dod   =exp(-deadtime*w*i);
Frek=gain*teller./nevner.*Dod;  % korrigert av SiS

%Amplitude og fase
A=abs(Frek);
vA=angle(Frek);

%Test for stor sprang
vA1=vA(1:length(vA)-1);
vA2=vA(2:length(vA));
DvA=vA2-vA1;

for ii=1:length(DvA)
  if DvA(ii) > 3
    vA(ii+1:length(vA))=vA(ii+1:length(vA))-2*pi;
  end
end


vA=vA*180/pi;

%Beregning av w180 og gain margin
if min(vA)>-180
  disp(['Ingen fase kryssover frekvens i det aktuelle frekvens' ...
	' intervallet'])
    w180=[];
    GM=[];
else
  for ii=1:length(vA)
    if vA(ii)>-180
      posisjon=ii;
    end
  end
  %Secant
  w180=w(posisjon)-(vA(posisjon)+180)* ... 
      (w(posisjon+1)-w(posisjon))/(vA(posisjon+1)-vA(posisjon));
  %Gain margin ved interpolasjon
  Lw180=A(posisjon)+(w180-w(posisjon))*(A(posisjon+1)-A(posisjon))...
      /(w(posisjon+1)-w(posisjon));
  GM=1/Lw180;
end



%Beregning av wc
if max(A)<1
  disp('Maks amplitude mindre enn ein.')
  wc=[];
  vAwc=[];
  PM=[];
else
  for ii=1:length(A)
    %Kryssing oven fra
    if A(ii)>1
      posisjon=ii;
    end
  end
  %Secant
  wc=w(posisjon)-(A(posisjon)-1)* ... 
      (w(posisjon+1)-w(posisjon))/(A(posisjon+1)-A(posisjon));
  %Fasen ved 
  vAwc=vA(posisjon)+(wc-w(posisjon))*(vA(posisjon+1)-vA(posisjon))...
      /(w(posisjon+1)-w(posisjon));
  PM=vAwc+180;
end

%Beregning av S og T
S=1./abs(1+Frek);
Ms=max(S);
T=abs(Frek./(1+Frek));
Mt=max(T);

if nargout==0
 figure
 title('Bode plot')
 subplot(211)
 loglog(w,A)
 hold
 loglog([w180,w180],[1,Lw180])
 loglog(w,ones(length(w)),':')
 ylabel('Amplitude')
 title('Bode plot')
 hold
 subplot(212)
 semilogx(w,vA)
 hold
 semilogx(w,-ones(length(w))*180,':')
 semilogx([wc,wc],[vAwc,-180])
 axis([min(w) max(w) max([min(vA),-360]) 0])
 ylabel('Fase forskyvning')
 xlabel('Frekvens')
 hold

 figure
 loglog(w,S,w,T,'.-')
 legend('S','T')
 xlabel('Frekvens')
 ylabel('Amplitude')
end


