function [ab,bb,cb,g,TL,TR] = balreal(a,b,c)
%BALREAL  Compute balanced state-space realization.
%
%   [SYSB,G] = BALREAL(SYS) computes a balanced state-space realization
%   of the LTI model SYS. For stable models, SYSB is an equivalent
%   realization for which the controllability and observability Gramians
%   are equal and diagonal, their diagonal entries forming the vector G
%   of Hankel singular values. This balances the input-to-state and
%   state-to-output energy transfers, and small entries of G indicate
%   states that can be removed with XELIM.
%
%   [SYSB,G] = BALREAL(SYS,Name1=Value1,Name2=Value2,...) specifies
%   additional options, see BALREALOPTIONS for details.
%
%   [SYSB,G,TL,TR] = BALREAL(SYS,...) also returns the balancing state
%   transformations TL,TR. For SYS with state x and matrices (A,B,C,D,E),
%   the balanced realization is (TL*A*TR,TL*B,C*TL,D,TL*E*TR) and the state
%   transformation is x = TR*x_b. Note that TL*E*TR=I for stable SYS.
%
%   Note:
%     * If SYS has unstable poles, its stable part is isolated, balanced,
%       and added back to the unstable part to form SYSB. The entries of G
%       corresponding to unstable modes are set to Inf.
%     * Use reducespec rather than BALREAL and XELIM for model reduction
%       purposes.
%
%   Example:
%     sys = zpk([-3 -.5],[-1e-6 -.4999 -10],1);
%     % Treat -1e-6 as unstable, notice small entry of G due to near
%     % pole/zero cancellation.
%     [sysb,g] = balreal(sys,ErrorBound='rel',Offset=.001)
%
%   See also BALREALOPTIONS, SSEQUIV, GRAM, XELIM, REDUCESPEC, SS.

% Old help
%warning(['This calling syntax for ' mfilename ' will not be supported in the future.'])
%BALREAL  Balanced state-space realization and model reduction.
%   [Ab,Bb,Cb] = BALREAL(A,B,C) returns a balanced state-space 
%   realization of the system (A,B,C).
%
%   [Ab,Bb,Cb,G,TL] = BALREAL(A,B,C) also returns a vector G containing
%   the diagonal of the gramian of the balanced realization, and 
%   matrix TL, the similarity transformation used to convert (A,B,C) 
%   to (Ab,Bb,Cb).  If the system (A,B,C) is normalized properly, 
%   small elements in gramian G indicate states that can be removed to
%   reduce the model to lower order.

%	J.N. Little 3-6-86
%	Revised 12-30-88
%       Alan J. Laub 10-30-94
%       P. Gahinet 6-27-96
%   Copyright 1986-2009 The MathWorks, Inc.

%       Reference:
%       [1] Laub, A.J., M.TL. Heath, C.C. Paige, and R.C. Ward,
%           ``Computation of System Balancing Transformations and Other
%           Applications of Simultaneous Diagonalization Algorithms,''
%           IEEE Trans. Automatic Control, AC-32(1987), 115--122.
%
%       The old balreal used an eigenvalue algorithm described in
%	 1) Moore, B., Principal Component Analysis in Linear Systems:
%	    Controllability, Observability, and Model Reduction, IEEE 
%	    Transactions on Automatic Control, 26-1, Feb. 1981.
%	 2) Laub, A., "Computation of Balancing Transformations", Proc. JACC
%	    Vol.1, paper FA8-E, 1980.

narginchk(3,3);
[sys,g,TR,TL] = balreal(ss(a,b,c,zeros(size(c,1),size(b,2))));
[ab,bb,cb] = ssdata(sys);
