function [Ar,Br,Cr,Tl,Tr] = balancedTruncation(A,B,C,r)
% BALANCEDTRUNCATION compute the balanced truncation reduced model
% for the linear system
%   \dot{x}(t) = Ax(t) + Bu(t)
%         y(t) = Cx(t)
% The algorithm computes matrices V,W\in\mathbb{R}^{n\times r}
% with W^TV = I_n and
%   Ar = W^TAV,     Br = W^TB,      Cr = CV
%
%
%
% WARNING: This code is only proof-of-concept and not tested or optimized
% for productional purposes!!!
%
%
%
%
% Input arguments:
%   A:      system matrix of size n times n
%   B:      input matrix of size n times m
%   C:      output matrix of size p times n
%   r:      the dimension of the reduced model
%
% Output arguments:
%   Ar:     reduced Ar
%   Br:     reduced Br 
%   Cr:     reduced Cr
%   Tl:     left projection matrix W
%   Tr:     right projection matrix V
%
% Authors:
%   Philipp Schulze, pschulze@math.tu-berlin.de
%   Benjamin Unger, unger@math.tu-berlin.de
%
% Version:
%   1.0 (December 22, 2017)
% 
PC          = lyapchol(A,B)' ;
QC          = lyapchol(A',C')' ;
[U,Sigma,V] = svd(PC'*QC) ;
U           = U(:,1:r) ;
V           = V(:,1:r) ;
Sigma       = Sigma(1:r,1:r) ;
Tl          = QC*V*diag(diag(Sigma).^(-0.5)) ;
Tr          = PC*U*diag(diag(Sigma).^(-0.5)) ;
Ar          = Tl'*A*Tr ;
Br          = Tl'*B ;
Cr          = C*Tr ;