interval_bode(NumVec,DenVec,clr)

function interval_bode(NumVec,DenVec,clr)
%interval_bode(NumVec,DenVec,clr)
%Plots the envelope of bode plot for interval systems.
%NumVec and DenVec are matrices
%with first row in each composed of lower limits and second row consisting
%of upper limits of the coefficients.
%
%Eg NumVec=[a4- a3- a2- a1- a0-;
%           a4+ a3+ a2+ a1+ a0+];
% if numerator is a4*s^4+a3s^3+a2s^2+a1*s+a0
%
% The argument clr indicates the colour with which the interval bode plot
% should be filled. It may be omitted if no filling is desired.
%
% Example :
% NumVec=[1 2 3 4 5 6;2 4 7 9 6 7];
% DenVec=[1 2 2 3 3 4 4;1 3 4 4 5 4 5];
%
% interval_bode(NumVec,DenVec,'b')

if (nargin < 2) || (nargin > 3)
    error('Function interval_bode should have 2 or 3 input arguments')
end
if ~isreal(NumVec) || ~isreal(DenVec)
    error('Numerator and Denominator Matrices should be real');
end

if (size(NumVec,1)~=2) || (size(DenVec,1)~=2)
    error('Numerator and Denominator Matrices should have exactly two rows');
end

if (any((NumVec(2,:)-NumVec(1,:)) < 0)) || (any((DenVec(2,:)-DenVec(1,:)) < 0))
    error('The elements in the second row should be greater than or equal to those in the first row of input vectors');
end

if nargin < 3
    clr='w';
else
if ~ismember(clr,'bgrcmykwBGRCMYKW')
    error('Invalid color argument')
end
end
Index=[1 1 -1 -1;
    1 -1 -1 1;
    -1 -1 1 1;
    -1 1 1 -1];
for i = 1:4
  
   N(i,:)= MakePoly(NumVec,Index(i,:));
   D(i,:)= MakePoly(DenVec,Index(i,:));
end
[m,p,w]=bode(tf(N(1,:),D(1,:)));
for i=1:4
    for j=1:4
        [m,p]=bode(tf(N(i,:),D(j,:)),w);
        l=4*i+j-4;
        m=reshape(m,[length(m) 1 1]);
        p=reshape(p,[length(p) 1 1]);
        mv(l,:)=20*log10(m)';
        pv(l,:)=p';
    end
end
magmax=max(mv);
magmin=min(mv);
phmax=max(pv);
phmin=min(pv);

subplot(2,1,1);
semilogx(w,magmax,w,magmin)
hold on
fill([w' w(end:-1:1)'],([magmax magmin(end:-1:1)]),clr)
ylabel('Magnitude (dB)')
title('Interval Bode Envelope')
set(gca,'XTickLabel',[])
%grid2 off x on y
subplot(2,1,2)
semilogx(w,phmax,w,phmin);
hold on
fill([w' w(end:-1:1)'],([phmax phmin(end:-1:1)]),clr)

ylabel('Phase (deg)')
xlabel('Frequency (rad/sec)')

function AnsPoly=MakePoly(IntervalPoly,Index)

Polylength=size(IntervalPoly,2);
PolyRep=ceil(Polylength/length(Index));
Index1=repmat(Index,[1 PolyRep]);
Index1=Index1(Polylength:-1:1);
AnsPoly=IntervalPoly(2,:);
I=find(Index1==-1);
AnsPoly(I)=IntervalPoly(1,I);
return