% ************************************************************************************************
% ********************** Program written by Stathis Skouras, July 2001 ***************************
% ************************************************************************************************

% running distillation model 'colamod_nonideal' for a ternary mixture
profile on

% Number of stages: 
NT=41;

% Normal boiling points of pure components (Reid et al, 1987)
TB1=337.7;    % = 63.9C, Component 1 (K) 
TB2=351.4;    % = 78.2C, Component 2 (K)
TB3=370.3;    % = 97.2C, Component 3 (K)

% States are : (NC-1)*NT mole fractions, NT liquid holdups, NT temperatures

% Initial conditions for states:
xA0=0.4*ones(1,NT);                           % Mole fraction of light component
xB0=0.4*ones(1,NT);                           % Mole fraction of intermediate component 
M0=0.4*ones(1,NT);                            % Liquid Holdup in top column section (kmol)  
T0 = xA0*TB1 + xB0*TB2 + (1-xA0-xB0)*TB3;     % Temperature at each stage (K)

% Defining options and mass matrix for DAE solver 
M=eye(3*NT,3*NT);                % 3*NT = no. of differential equations (compositions, liq.holdup)
M(3*NT+1:4*NT,3*NT+1:4*NT)=0;    % NT   = no. of algebraic equations (pressure -> temperature)
options.Mass=M; 

% Running column
[t,x]=ode15s(@cola_lv_nonideal,[0 20000],[xA0,xB0,M0,T0],options);

% Save results
lengthx=size(x); 

% Vector 'Xss' returns steady-state molar compositions of component A and B at each stage,
% molar liquid holdups at each stage and temperature at each stage 
Xss=x(lengthx(1),:)'

profile off
profile report
