% script for Question 8/Leuven 2005Gbal,sig]
N=100;
A=zeros(2*N);
B=zeros(2*N,1);
C=zeros(1,2*N);
for k=1:N,
    k2=2*k;
    k1=k2-1;
    A(k1:k2,k1:k2)=[0 1;-k/10 -2];
    B(k1:k2)=[0;1];
    C(k1:k2)=[1 0];
end
G=ss(A,B,C,0);
Gnom=pck(A,B,C,0);
[Gbal,sig]=sysbal(Gnom);
Gred=hankmr(Gbal,sig,10,'d');
[Ar,Br,Cr,Dr]=unpck(Gred);
Gr=ss(Ar,Br,Cr,Dr);
fprintf('\nq8leuven: lower bound %d, actual %d\n',sig(11),norm(G-Gr,Inf));
ww=linspace(0,10,300);
g=squeeze(freqresp(G,ww));
gr=squeeze(freqresp(Gr,ww));
close(gcf);
subplot(2,1,1); plot(ww,real(g),ww,real(gr));grid
subplot(2,1,2); plot(ww,imag(g),ww,imag(gr));grid