clear all;
load('cola_init.mat');
global T_l T_r match_full salva
[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;
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));
for i=73:length(g)
    pos_keep_g(i) = 0;
end
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;

match_full = x_SS(end,:)';
x_reduced = fsolve(@fun_sol,x_SS(end,:)');
z_init = T_l*x_reduced;
tic 
[t_DYN_red,z_DYN_red] = ode15s('cola_lv_F1_red',[t_SS(end) t_SS(end)+1200],z_init);
t_72 = toc;
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);




pos_keep_g = ones(length(g));
for i=80: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);

match_full = x_SS(end,:)';
x_reduced = fsolve(@fun_sol,x_SS(end,:)');
z_init = [];
z_init = T_l*x_reduced;
tic
[t_DYN_red_80,z_DYN_red_80] = ode15s('cola_lv_F1_red',[t_SS(end) t_SS(end)+1200],z_init);
t_80 = toc;
for i=1:length(t_DYN_red_80)
    x_DYN_red_80(i,:) = T_r*z_DYN_red_80(i,:)';
end

plot_t_RED_80 = [t_DYN_red_80];
plot_x_RED_80 = [x_DYN_red_80];
plot_xB_RED_80 = plot_x_RED_80(:,1); 
plot_yD_RED_80 = plot_x_RED_80(:,41);






figure(1)
plot(plot_t_FULL',plot_yD_FULL');
hold on
plot(plot_t_RED',plot_yD_RED','r');
plot(plot_t_RED_80',plot_yD_RED_80','*r');
hold off
title('concentration vapour top product');
legend('FULL, 82 states','REDUCED, 72 states','REDUCED, 79 states');
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');
plot(plot_t_RED_80',plot_xB_RED_80','*r');
hold off
title('concentration liquid bottom product');
legend('FULL, 82 states','REDUCED, 72 states','REDUCED, 79 states');
xlabel('time [min]');
ylabel('x_B [-]');

%%
