%
%   This script examplifies Kalman decomposition in Matsys Gk
%
%   Given system
A = [ 2 4 3; 4 2 3; 4 4 4];
B = [1 -1 0; 1 1 -2]';
C = [2 2 -3];

%   Reachability
gamma = ctrb(A,B)

%   Observability
omega = obsv(A,C)

%   Basis for reachable and unobservable
e1 = [ 1 -1 0 ]';
%   Basis for reachable and observable 
e2 = [ 1 1 -2 ]';
%   Basis for unreachable and observable
e3 = [ 0 3 2 ]';

%   Change of variables
T = [ e1 e2 e3 ];

%   Transformed system
Atilde = inv(T)*A*T
Btilde = inv(T)*B
Ctilde = C*T

% Check of minimal realization
sys = ss(A,B,C,0);
zpk(sys)
sysmin = ss(Atilde(2,2),Btilde(2,:),Ctilde(:,2),0);
zpk(sysmin)

% >> kalmandecomposition
% 
% gamma =
% 
%      1     1    -2     0     4     0
%     -1     1     2     0    -4     0
%      0    -2     0     0     0     0
% 
% 
% omega =
% 
%      2     2    -3
%      0     0     0
%      0     0     0
% 
% 
% Atilde =
% 
%    -2.0000   -0.0000   18.0000
%          0         0   -0.0000
%          0         0   10.0000
% 
% Btilde =
% 
%      1     0
%      0     1
%      0     0
% 
% Ctilde =
% 
%      0    10     0
% 
%  
% Zero/pole/gain from input 1 to output:
% 1.7764e-15 (s-10)
% -----------------
%  s (s-10) (s+2)
%  
% Zero/pole/gain from input 2 to output:
% 10 (s-10) (s+2)
% ---------------
% s (s-10) (s+2)
%  
%  
% Zero/pole/gain from input 1 to output:
% 0
%  
% Zero/pole/gain from input 2 to output:
% 10
% --
% s