load build
% load CD
% load beam
% load iss   % slow decay

 A = full(A);
 B = full(B);
 C = full(C);
 n = length(A)

 [Us,A] = schur(A);
 B = Us'*B;
 C = C*Us; 


 U = lyapUR(A,B);

 id = [n:-1:1];
 L = lyapUR(A(id,id)',C(id)');
 L = L(id,id);

 P = U*U';
 Q = L*L';
 fprintf('norm(A*P+P*A''+B*B'') = %8.5e\n', norm(A*P+P*A'+B*B'))
 fprintf('norm(A''*Q+Q*A+C''*C) = %8.5e\n', norm(A'*Q+Q*A+C'*C))

 figure(1), clf
 svP = svd(P);
 svQ = svd(Q);
 Hsv = svd(P*Q);
 semilogy(svP,'b.','markersize',12), hold on
 semilogy(svQ,'r.','markersize',12)
 semilogy(Hsv,'k.','markersize',12)
 set(gca,'fontsize',18)
 legend('singular values of P','singular values of Q', 'Hankel singular values')
 xlabel('$k$','fontsize',16,'interpreter','latex')
 ylabel('$\sigma_k$','fontsize',16,'interpreter','latex')
 

 [Z,Sig,Y] = svd(U'*L);
 
 sig = diag(Sig);
 W = L*Y*diag(1./sqrt(sig));
 V = U*Z*diag(1./sqrt(sig));
 figure(2), clf
 modn=Mysigma_log(A,B,C,D,-2,7,'Frequency Response','k-'); hold on
 xlim([1e-2 1e7])
 pause

 clear modk
 set(gca,'fontsize',18)
 for k=1:n
     Ak = W(:,1:k)'*A*V(:,1:k);    
     Bk = W(:,1:k)'*B;
     Ck = C*V(:,1:k);    
     if exist('modk'), delete(modk), end
     modk = Mysigma_log(Ak,Bk,Ck,D,-2,7,'Frequency Response','r-');
     legend([modn modk],sprintf('full, n=%d',n), sprintf('k=%d',k),3)
     pause
 end