function K=ps54(T,n,d)
% function K=ps54(T,n,d)
%
% T  -  delay (same as tau)
% n  -  initial approximation degree (n<4 means using L=1/s)
% d  -  small parameter

if nargin<1, T=1; end
if nargin<2, n=1; end
if nargin<4, d=1e-3; end

s=tf('s');
P=(s-1)/(s*(s-2));
m=min(3,max(0,n-3));
[A,B,C,D]=ssdata(ss((1-s*(0.5*T/n))/(1+s*(0.5*T/n)))^n);
sys=pck(A,B,C/A,0);
[sys,sig]=sysbal(sys);
[Ar,Br,Cr,Dr]=unpck(hankmr(sys,sig,m,'d'));
L0=ss(Ar,Br,Cr,Dr);
%L0=hankelmr(ss(A,B,C*inv(A),0),m);
L=1/s+L0;
ww=logspace(-2,2,10000);
sss=1i*ww(:);
mm=abs(squeeze(freqresp(L0,sss))-(exp((-T)*sss)-1)./sss);
r=max(mm);
     
assignin('base','d',d)         % export variables
assignin('base','P',P)
assignin('base','L',L)

p=linmod('ps54des');           % extract open loop
p=ss(p.a,p.b,p.c,p.d);
K=hinfsyn(p,1,1);              % optimize controller

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