%Illustration of admissible process uncertainty
%The nominal loop transfer function is plotted
%together with the uncertainty circles
%kja 040826

aminit;					% initialize MATALB

figure(1); clf;
clear i;

w=logspace(-5, 1, 1000); s=i*w;
b =1.3785; a=-0.0142; k=0.7151; ki=0.1814;
P=b./(s-a);
L=P.*(k+ki./s);
T=L./(1+L);
T0=T;
[Mt,iii]=max(abs(T));
wmt=w(iii);
Mt,wmt
absP=abs(P(iii));
absPoT=abs(P(iii)/T(iii));
absP,absPoT
[xxx,jjj]=min(abs(w-5));
w2=w(jjj);
Tw2=abs(T(jjj));
w2,Tw2
1/Tw2

% Plot the complementary sensitivity function and the max delta P
subplot(221); 
loglog([0.1 10], [1 1], 'k-', 'LineWidth', AM_ref_linewidth); hold on;
loglog([1 1], [0.1 10], 'k:', 'LineWidth', AM_ref_linewidth);
%pl=loglog([min(w), max(w)], [1 1], ...
%  'r--', w, 1./abs(T), 'g-.', ...
%  w, abs(P) ./ abs(T), 'b-');
pl=loglog(w, 1./abs(T), '-', w, abs(P) ./ abs(T), '--'); hold on;
set(pl,'Linewidth', AM_bode_linewidth);
amaxis([0.1 10 0.5 10]);

set(gca, 'XMinorGrid', 'off', 'YMinorGrid', 'off');
xlabel('freq');
ylabel('unc');

% Add some labels for what we are looking at
amtext(0.27, 6, 'P/T');
amtext(3, 6, '1/T');

% Plot the Nyquist curve
subplot(222);
ampolezero(tf([1], [1]), [-50 225 -125 125]);
np=plot(real(P),imag(P),'b'); set(np,'Linewidth', AM_data_linewidth)
hold on

% set up the color
color = 'r';
if (AMPRINT_FLAG==2) color = 'k'; end;

% Now generate the uncertainty circles
%wp = [10^-12 0.1 0.18 0.3 0.65 1 2];
wp = [10^-12 abs(a) 0.05];
th0 = pi/4;
th = -pi/3;
for ind = 1:length(wp)
  s=i*wp(ind); P=b./(s-a);
  L=P.*(k+ki./s); T=L./(1+L); S=1./(1+L);
  cx=real(P); cy=imag(P); r=abs(P)/abs(T);
  circplot(cx,cy,r,'r--', 'LineWidth', AM_ref_linewidth);
  plot(cx, cy, 'ro', 'MarkerSize', AM_data_markersize);
  arrow([cx, cy], [cx, cy] + r*[cos(th*ind + th0), sin(th*ind + th0)], ...
    AM_normal_arrowsize, 'FaceColor', color, 'EdgeColor', color);
end

amaxis([-50 200 -125 125]);
% text(2.5,-0.2,'x');
% text(0.5,2,'y');
xlabel('x'); ylabel('y');
axis square; 
hold off;

amprint('cruise-robustness.eps');

figure(2); clf;
P = tf([b], [1 -a]);
C = tf([k ki], [1 0]);
L = P*C;
bodemag(P, C, L, P*(1 + L)/L)
legend('P', 'C', 'L' ,'P/T')