function [K,S]=ps42(a,b,d)
% function [K,S]=ps42(a,b,d)
%
% 
if nargin<1, a=1; end          % default a
if nargin<2, b=0.43*a; end     % default requested bandwidth
if nargin<3, d=0.001/a; end    % default small weight

s=tf('s');                     % useful constant
P0=(s-a)/(s*(s+2*a));          % the plant
[A,B,C,D]=butter(3,b,'s');     % Butterworth filter
W=ss(A,B,C,D);
W=(10/abs(squeeze(freqresp(W,1i*b))))*W;
%ww=linspace(0,2*b,500);
%gg=abs(squeeze(freqresp(W,1i*ww)));
%close(gcf);plot(ww,abs(gg));grid

assignin('base','d',d)         % export variables
assignin('base','P0',P0)
assignin('base','W',W)

p=linmod('ps42des');           % extract open loop
p=ss(p.a,p.b,p.c,p.d);
[K,G]=hinfsyn(p,1,1);          % optimize controller
fprintf('success flag for w0=%f:  %f<1\n',b,norm(G,Inf))

assignin('base','K',K)         % export controller
S=linmod('ps42test');          % extract closed loop
fprintf('max(real(eig(A)))=%f\n',max(real(eig(S.a))))
S=ss(S.a,S.b,S.c,S.d);
fprintf('norm(S,Inf)=%f\n',norm(S,Inf))
ww=logspace(-2,5,10000);
Sw=abs(squeeze(freqresp(S,1i*ww)));
Kw=abs(squeeze(freqresp(K,1i*ww)));
k=find(Sw>0.1,1);              % find actual bandwidth
w1=ww(k-1);
fprintf('actual bandwidth:  %f (vs. %f)\n',w1,b)
close(gcf);
subplot(2,1,1);loglog(ww,Sw,[w1 w1],[0.01 0.1], ...
    [ww(1) ww(end)],[0.1 0.1],[ww(1) ww(end)],[20 20]);grid
xlabel('w (rad/sec)')
ylabel('|S(jw)|')
subplot(2,1,2);loglog(ww,Kw);grid
xlabel('w (rad/sec)')
ylabel('|K(jw)|')