function H=hsvd_6242(G,m)
% function H=hsvd_6242(G,m)
%
% Hankel svd for CT stable system G
%   with one argument:  H is the ordered vector of Hankel singular values
%   with two arguments: H is the m-th order btr reduced system

if nargin<1,   % a test example
    G=ss(diag(-(1:30)),ones(30,1),ones(1,30),0);
    m=3;
end

[A,B,C,d]=ssdata(G);
Wc=lyap(A,B*B');
Wo=lyap(A',C'*C);
[V,D]=eig(Wc*Wo);
V=real(V);
[W,I]=sort(-sqrt(abs(diag(D))));
V=V(:,I);
W=-W;
if nargin==1,
    H=W;
else
    V=V(:,1:m);
    U=(V'*Wo*V)\(V'*Wo);
    H=ss(U*A*V,U*B,C*V,d);
end

if nargin~=1, 
    w=linspace(0,100,1000);
    g=squeeze(freqresp(G,w));
    g1=squeeze(freqresp(H,w));
    close(gcf)
    subplot(2,1,1);plot(w,real(g),w,real(g1)); grid
    subplot(2,1,2);plot(w,imag(g),w,imag(g1)); grid
end