clear all
A = [-2.0 0.0 0.0; 1.0 -1.1 0.0; 0.0 0.1 -1.0];
B = [2.0; 0.0; 0.0];
C = [0.0 0.0 1.0];
D = 0.0;

sys = ss(A,B,C,D);

[sysb,g,T,Ti] = balreal(sys);

pos_keep_g = (g> 1e-36);
T_l = T;
T_r = Ti;


for i=length(g):-1:1
    if (pos_keep_g(i)==0)
        
        T_l(i,:) = [];
        
        T_r(:,i) = [];
    end;
end

new_dym = size(T_l,1);
% x = lyap(a,b,c) solves
% ax + xb + c = 0

% A*P + P*A^T + B*B^T = 0
P = lyap(A,A',B*(B'));

% A^T * Q + Q*A + C^T * C = 0
Q = lyap(A',A,(C')*C);

% r = lyapchol(a,b) solves
% a*x + x*a' + b*b' = 0 with x =r'r

% A*P + P*A^T + B*B^T = 0
%R_P = lyapchol(A,B);
%P = R_P'*R_P;

% A^T * Q + Q*A + C^T * C = 0
%R_Q = lyapchol(A',C');
%Q = R_Q'*R_Q;

%[t_DYN_red,z_DYN_red] = ode15i('cola_lv_F1_red_i',[t_SS(end) t_SS(end)+1200],T_l*(x_SS(end,:)'),zeros(new_dym,1));
[t_DYN_red,z_DYN_red] = ode15s('cola_lv_F1_red',[t_SS(end) t_SS(end)+1200],T_l*(x_SS(end,:))');

for i=1:length(t_DYN_red)
    x_DYN_red(i,:) = T_r*z_DYN_red(i,:)';
end

plot_t_RED = [t_DYN_red];
plot_x_RED = [x_DYN_red];
plot_xB_RED = plot_x_RED(:,1); 
plot_yD_RED = plot_x_RED(:,41);

figure(1)
plot(plot_t_FULL',plot_yD_FULL');
hold on
plot(plot_t_RED',plot_yD_RED','r');
hold off
title('concentration vapour top product');
legend('FULL','REDUCED');
xlabel('time [min]');
ylabel('y_D [-]');

figure(2)
plot(plot_t_FULL',plot_xB_FULL');
hold on
plot(plot_t_RED',plot_xB_RED','r');
hold off
title('concentration liquid bottom product');
legend('FULL','REDUCED');
xlabel('time [min]');
ylabel('x_B [-]');