function [Ak,Bk,Ck,s] = balS(A,B,C,tol);
%
%
%     Usage:  [Ak,Bk,Ck,s] = balS(A,B,C,tol);
%
%     Produces a truncated balanced realization of the LTI state 
%     space system (A,B,C).  
%
%     A is an n by n matrix 
%     with all eigenvalues in the open left halfplane.
%
%     B is an n by m matrix (m << n)
%     C is an p by n matrix (p << n)
%
%     This routine calls lyapSL and lyapSU which compute
%     Cholesky factors of system Grammians P and Q 
%     This is done in the real Schur basis.  The only complex 
%     arithmetic used is to compute eigensystem of 2 by 2 matrices
%     when there is a 2 by 2 block on the diagonal of the quasi upper 
%     triangular matrix R, the real Schur form of A.
%
%     D.C. Sorensen
%     sorensen@rice.edu
%     8 Jan 00
%
%--------------------------------------------------------------
%
%
%  Compute the real Schur decomposition of A.
%

    sys  = ss(A,B,C,0);
    [Q,R] = schur(A);
%
%  initialize
%
    C = C*Q; B = Q'*B;
    [n,n] = size(A);
%
%  Solve for L in the Real Schur basis.  
%
    [L] = lyapSL(R,C);
    M = R'*L;
    Lnorm = norm(M*L' + L*M' + C'*C,1)/(norm(L,1)^2)

%
%  Solve for U in the Real Schur basis.  
%
    [U] = lyapSU(R,B);
    M = R*U;
    Unorm = norm(M*U' + U*M' + B*B',1)/(norm(U,1)^2)


   [W,S,Y] = svd(U'*L); 


%  A1 = (Y'*(L'*R*U)*W);
% 
%  B = Y'*(L'*B);
%  C = (C*U)*W;
%
%  M1 = A1*S + S*A1' + B*B';
%  resU = norm(M1,1)/(norm(B,1)^2)
%
%  M2 = A1'*S + S*A1 + C'*C;
%  resL = norm(M2,1)/(norm(C,1)^2)

%
%  Determine first k such that  s(k+1) < tol*s(1)
%
   s = diag(S);
   k = 1;
   while(k < n & s(k) > tol*s(1)), k = k+1; end

%
%  Truncate and balance the system
%
   Sk = sqrt(S(1:k,1:k));
   Ak = Sk\(L*Y(:,1:k))'*R*(U*W(:,1:k))/Sk;
   Bk =Sk\(L*Y(:,1:k))'*B;
   Ck= C*(U*W(:,1:k))/Sk;

   Sk = diag(s(1:k));

%  Bnorm = norm(Ak*Sk + Sk*Ak' + Bk*Bk',1)/(norm(Bk,1)^2)
%  Cnorm = norm(Ak'*Sk + Sk*Ak + Ck'*Ck,1)/(norm(Bk,1)^2)
%
%  display graphics of sigma plots
%

   sysk = ss(Ak,Bk,Ck,0);
         
   subplot(2,1,1), sigma(sys,'b',sysk,'r');         
 

   sysE = sys - sysk;
   subplot(2,1,2), sigma(sysE)


   
