%% Main program
clear all; clc; close all

parameters_c

% Inputs and disturbances:
LT = 226.8/60;          % Reflux [kmol/min]
VB = 309.4/60;          % Vapour boilup [kmol/min]
F  = 350/60;            % Feedrate [kmol/min]
qF = 1.064;             % Feed liquid fraction
%     MetOH     TG  DG MG GL BD
zF = [333/350 17/350 0 0 0 0]; % Feed molar compositions
u1 = [LT VB];
d  = [F qF zF];

% Initial values:
load steadystate_5 % y0
% load steadystate_high_350
T0 = (91.27+273.15)*ones(NT,1);
x0 = [y0; T0];

%% Steady-state calculations
options = optimset('algorithm','interior-point','Display','iter','Diagnostics','on');
tic
x  = fmincon(@cost_c,x0,[],[],[],[],[],[],@reactivedist_c,options,u1,d);
toc
f1 = output_c(x,u1,d);
print_results_c(f1)

%% Calculating the conversion
xB = f1{1,6};
xTGB = xB(1,2);
FTG = xTGB*f1{1,4};
F0TG = F*zF(1,2);
con = (F0TG - FTG)/F0TG
BBD = f1{1,4}*(1-sum(xB));
yield = BBD/(3*F0TG)

x0 = x;
%% Optimization calculations
options = optimset('algorithm','interior-point','Display','iter','Diagnostics','on','maxfunevals',150000, 'maxiter',5000);
lb = zeros(length([x0;u1']),1);
tic
y  = fmincon(@cost_c,[x0;u1'],[],[],[],[],lb,[],@reactivedist_c,options,d);
toc
z0 = y(1:length(x0)); u2 = y(length(x0)+1:end);
f2 = output_c(z0,u2,d);
print_results_c(f2);

LT = f2{1};       % Reflux
VB = f2{2};       % Boilup
u1 = [LT VB];

%% Dynamic calculations
% Using DAE system and mass matrix
% F  = 350;       % Introducing disturbance in feed
% d  = [F qF zF];     % Disturbances
% 
% % The mass matrix:
% temporary = ones(NT*NC,1);
% temporary1 = zeros(NT,1);
% m_matrix = diag([temporary;temporary1]);
% options = odeset('mass', m_matrix);
% 
% % Time:
% t0 = 0;
% tfin = 10;
% 
% % Solving the system:
% tic
% [t,zf] = ode15s(@reactivedist_c_DAE,[t0 tfin],z0,options);
% toc
% f3 = output_dynamic(zf(end,:),u1,d);