function ps83(g,m,n,dd)
% function ps83(g,m,n,dd)
%
% g  - desired closed loop L2 gain bound
% m  - desired controller order
% n  - order of the initial approximation for (exp(-s)-exp(-1))/(s-1)
% dd - vector of scaling factors to use

if nargin<1, g=100; end
if nargin<2, m=2; end
if nargin<3, n=100; end
if nargin<4, dd=logspace(-1,1,10)'; end
dd=dd(:);
nd=length(dd);
s=tf('s');
[A,B,C,D]=ssdata(ss((1-s*(0.5/n))/(1+s*(0.5/n)))^n);
sys=pck(A,B,C/(A-eye(n)),0);
[sys,sig]=sysbal(sys);
[Ar,Br,Cr,Dr]=unpck(hankmr(sys,sig,m-1,'d'));
L0=ss(Ar,Br,Cr,Dr);    % approximation for (exp(-s)-exp(-1))/(s-1)
%L0=hankelmr(ss(A,B,C/(A-eye(n)),0),m);
P0=exp(-1)/(s-1)+L0;
ww=logspace(-2,2,10000);
sss=1i*ww(:);
mm=abs(squeeze(freqresp(L0,sss))-(exp(-sss)-exp(-1))./(sss-1));
r=max(mm);
%close(gcf);semilogx(ww,mm,[ww(1) ww(end)],[r r]);grid;pause
fprintf('approximation error: %f\n',r)

h=sqrt(g*r);   
assignin('base','P0',P0)
assignin('base','h',h)
gg=zeros(nd,1);
for i=1:length(dd),
    d=dd(i);
    assignin('base','d',d)         % export variables
    p=linmod('ps83des');           % extract open loop
    p=ss(p.a,p.b,p.c,p.d);
    [K,~,GAM]=hinfsyn(p,1,1);      % optimize controller
    gg(i)=GAM;
    fprintf('.')
end
Af=ssdata(K);
fprintf('order: %d\n',size(Af,1))
close(gcf); plot(dd,gg/g); grid
return

assignin('base','K',K)         % export controller
p=linmod('ps54test');          % extract closed loop
S=ss(p.a,p.b,p.c,p.d);
fprintf('[T=%f]  stability: %f<0,  small gain: %f<%f\n', ...
    T,max(real(eig(p.a))),norm(S,Inf),1/r)
close(gcf);plot(ww,mm,[ww(1),ww(end)],[r r]);grid