% ranges

ms=.9:.1:1.1;
bs=.1:.05:.2;
ks=2:.5:3;

% nominal

m0=1;
b0=.15;
k0=2.5;

P0=tf(1,[m0,b0,k0]);

% Plot all TFs

sys=cell(length(ms)*length(bs)*length(ks),1);
i=1;
for m=ms
  for b=bs
    for k=ks    
      sys{i}=tf(1,[m,b,k]);
      i=i+1;      
    end
  end
end

figure(1);clf
bode(sys{:},{.1,20})
grid on

print -depsc '../tf_range.eps'

title('Process Bode plot')
disp('press <return> to continue');pause

% Plot all additive bounds

sys=[];
i=1;
for m=ms
  for b=bs
    for k=ks    
      sys=[sys;tf(1,[m,b,k])-P0];
      i=i+1;      
    end
  end
end

[m,p,w]=bode(sys,{.1,20});
m=20*log10(reshape(m,size(m,1),size(m,3)));
figure(2);clf
semilogx(w,m);
axis([.1,20,-100,20])
xlabel('Magnitude (dB)')
ylabel('Frequency (rad/sec)')
grid on
ladd=[.01,.6,1.3,1.9,4,20;
     -14,-14,19,19,-32,-66];
line(ladd(1,:),ladd(2,:),'marker','o','linewidth',3,'linestyle','--')

print -depsc '../add_bound.eps'

title('Additive uncertainty')
disp('press <return> to continue');pause


% Plot all multiplicative bounds

sys=[];
i=1;
for m=ms
  for b=bs
    for k=ks    
      sys=[sys;tf(1,[m,b,k])/P0-1];
      i=i+1;      
    end
  end
end

[m,p,w]=bode(sys,{.1,20});
m=20*log10(reshape(m,size(m,1),size(m,3)));
figure(3);clf
semilogx(w,m);
axis([.1,20,-100,20])
xlabel('Magnitude (dB)')
ylabel('Frequency (rad/sec)')
grid on
lmult=[.01,.6,1.3,1.9,4,100;
      -6,-6,19,19,-12,-12];
line(lmult(1,:),lmult(2,:),'marker','o','linewidth',3,'linestyle','--')

print -depsc '../mult_bound.eps'

title('Multiplicative uncertainty')
disp('press <return> to continue');pause

%% Controller 1

% Integral
C=zpk([],0,.25);
Ceps=zpk([],-1e-2,.25);

% rltool(P0,C)

% PID
C=zpk([-1,-2],0,5);
Ceps=zpk([-1,-2],-1e-2,5);

%% Nominal nyquist

figure(4);clf
nyquist(Ceps*P0);


print -depsc '../nyquist_big.eps'

title('Nyquist plot for 1st controller')
disp('press <return> to continue');pause

axis([-20,30,-20,20])

print -depsc '../nyquist_small.eps'


%% Nominal sensitivity

S0=feedback(1,C*P0);

figure(5);clf
bode(1-S0)
grid on

ax=get(gcf,'children');
axes(ax(2));
line(lmult(1,:),-lmult(2,:),'marker','o','linewidth',3,'linestyle','--')
ax=axis;
ax(3)=min([ax(3),-lmult(2,:)]);
ax(4)=max([ax(4),-lmult(2,:)]);
axis(ax);

print -depsc '../small_gain1.eps'

title('Small-gain condition for 1st controller')
disp('press <return> to continue');pause

%% Controller 2

C   =zpk([-5,-5],[0    ,-.5,-.5],.005);
Ceps=zpk([-5,-5],[-1e-2,-.5,-.5],.005);

%rltool(P0,C)

S0=feedback(1,C*P0);

figure(6);clf
bode(1-S0)
grid on

ax=get(gcf,'children');
axes(ax(2));
line(lmult(1,:),-lmult(2,:),'marker','o','linewidth',3,'linestyle','--')
ax=axis;
ax(3)=max(min([ax(3),-lmult(2,:)]),-60);
ax(4)=max([ax(4),-lmult(2,:)]);
axis(ax);

print -depsc '../small_gain2.eps'

title('Small-gain condition for 2nd controller')
disp('press <return> to continue');pause