clear all;
load('cola_init.mat');
global T_l T_r T Ti z_init_2
[t_SS,x_SS] = ode15s('cola_lv',[0 100],Xinit);

tic
[t_DYN,x_DYN] = ode15s('cola_lv_F1',[t_SS(end) t_SS(end)+1200],x_SS(end,:)');
t_FULL = toc;
sprintf('FULL MODEL, CPU time: %f', t_FULL)

plot_t_FULL = [t_SS; t_DYN];
plot_x_FULL = [x_SS; x_DYN];
plot_xB_FULL = plot_x_FULL(:,1); 
plot_yD_FULL = plot_x_FULL(:,41);

t_DYN = [];
x_DYN = [];

Ls=2.70629; Vs=3.20629; Fs=1.0; zFs=0.5;
[A,B,C,D]=cola_linearize('cola_lv_lin_ivan',x_SS(end,:),[Ls Vs Fs zFs]);

sys = ss(A,B,C,D);
[sysb,g,T,Ti] = balreal(sys);

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

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);

z_init = [];
z_init = T*x_SS(end,:)';
z_init_1 = [];
z_init_1 = z_init(1:new_dym,1);
z_init_2 = [];
z_init_2 = z_init(new_dym+1:end,1);

t_DYN_red = [];
z_DYN_red_1 = [];
x_DYN_red = [];
tic
[t_DYN_red,z_DYN_red_1] = ode15s('cola_lv_F1_red',[t_SS(end) t_SS(end)+1200],z_init_1);
t_RED = toc;
sprintf('RED MODEL caso 1, number state:%d CPU time: %f', new_dym ,t_RED)

z_DYN_red = [];
for i=1:length(t_DYN_red)
    z_t = [];
    z_t = [z_DYN_red_1(i,:), z_init_2'];
    z_DYN_red = [z_DYN_red; z_t];
    x_DYN_red(i,:) = Ti*z_DYN_red(i,:)';
end

plot_t_RED_caso_1 = [t_DYN_red];
plot_x_RED_caso_1 = [x_DYN_red];
plot_xB_RED_caso_1 = plot_x_RED_caso_1(:,1); 
plot_yD_RED_caso_1 = plot_x_RED_caso_1(:,41);



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

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);

z_init = [];
z_init = T*x_SS(end,:)';
z_init_1 = [];
z_init_1 = z_init(1:new_dym,1);
z_init_2 = [];
z_init_2 = z_init(new_dym+1:end,1);

t_DYN_red = [];
z_DYN_red_1 = [];
x_DYN_red = [];
tic
[t_DYN_red,z_DYN_red_1] = ode15s('cola_lv_F1_red',[t_SS(end) t_SS(end)+1200],z_init_1);
t_RED = toc;
sprintf('RED MODEL caso 2, number state:%d CPU time: %f', new_dym ,t_RED)

z_DYN_red = [];
for i=1:length(t_DYN_red)
    z_t = [];
    z_t = [z_DYN_red_1(i,:), z_init_2'];
    z_DYN_red = [z_DYN_red; z_t];
    x_DYN_red(i,:) = Ti*z_DYN_red(i,:)';
end

plot_t_RED_caso_2 = [t_DYN_red];
plot_x_RED_caso_2 = [x_DYN_red];
plot_xB_RED_caso_2 = plot_x_RED_caso_2(:,1); 
plot_yD_RED_caso_2 = plot_x_RED_caso_2(:,41);


figure(1)
plot(plot_t_FULL',plot_yD_FULL');
hold on
plot(plot_t_RED_caso_1',plot_yD_RED_caso_1','r');
plot(plot_t_RED_caso_2',plot_yD_RED_caso_2','*r');
hold off
title('concentration vapour top product');
legend('FULL','REDUCED caso 1','REDUCED caso 2');
xlabel('time [min]');
ylabel('y_D [-]');

figure(2)
plot(plot_t_FULL',plot_xB_FULL');
hold on
plot(plot_t_RED_caso_1',plot_xB_RED_caso_1','r');
plot(plot_t_RED_caso_2',plot_xB_RED_caso_2','*r');
hold off
title('concentration liquid bottom product');
legend('FULL','REDUCED caso 1','REDUCED caso 2');
xlabel('time [min]');
ylabel('x_B [-]');

%%
% 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;
