
% Initializes & runs the model

clear all; clc;

global coldata dyndata hydraulic termodata U NI tend

%% COLUMN DATA %%
coldata.Nt = [7 10 4 4 4 5 4];

% Initial holdup in reboiler and condenser
coldata.M0 = [0.5 0.5 1 0.5 0.5 1 0.5];
hydraulic.k1 = 1/0.063;

%% THERMODYNAMIC DATA %%

% Universal gas constant [cal/molK]
termodata.R = 1.98721;
% Antoine coefficients
termodata.A = [
    8.08097  1582.271  239.726
    8.11220  1592.864  226.184
    8.37895  1788.020  227.438
    7.36366  1305.198  173.427
    ];
% Molar volumes of pure components
termodata.V = [40.73  58.68  75.14  91.97];
% Wilson parameters
termodata.W = [
    0           507.3868    71.4631     -121.3194
    -452.5062   0           102.6706    -8.3924
    154.6476    -40.8745    0           -316.7718
    9543.2503   122.8635    547.4636    0
    ];
% Pressure (atm)
termodata.P = 1.0;
%%
indmat

U=[]; NI=9;


load Nominal

dyndata.x_nom=x_nom;
dyndata.u_nom=u_nom;


options = odeset('RelTol',1e-13,'AbsTol',1e-10);
tend=20000;
[t,x1]=ode15s(@Cascade,[0 tend],[zeros(NI,1);x_nom],options);

x=x1(end,NI+1:end)';
u=U(end,2:end);
dxdt=Kaibel_dxdt(0,[u';x]);
NORM_dxdt=norm(dxdt)
u=U(end,2:end)';

CheckSum
[[x(im_xt());1-sum(x(im_xt()))] [x(im_xS1());1-sum(x(im_xS1()))] [x(im_xS2());1-sum(x(im_xS2()))] [x(im_xr());1-sum(x(im_xr()))]]


