% running distillation model 'multivessel' for a ternary mixture
% Number of stages in the top section: 
NT1=25;
%Number of stages in the bottom section (without the reboiler)
NT2=25; 
% Total number of 'stages'  
NT=NT1+NT2+3;                             % We include condenser, middle vessel and reboiler 
    
% Initial conditions
zFA=0.6*ones(1,NT);                           % Initial composition of light component
zFB=0.2*ones(1,NT);                           % Initial composition of intermediate component
Md=0.5;                                     % Initial holdup in condenser (kmol)
M1=1/500*ones(1,NT1);                         % Initial holdup in top column section (kmol)  
Mm=2.0;                                     % Initial holdup in middle vessel (kmol)
M2=1/500*ones(1,NT2);                         % Initial holdup in bottom column section (kmol) 
Mb=2.5;                                      % Initial holdup in reboiler (kmol) 
T0=(75+273.15)*ones(1,NT);                   % Initial temperature in each stage (K)

% Defining options and mass matrix for ODE solver 
M=eye(3*NT,3*NT);
M(3*NT+1:4*NT,3*NT+1:4*NT)=0;
options.Mass=M; 

% Running column
[t,x]=ode15s(@input_multivessel,[0 1],[zFA,zFB,Md,M1,Mm,M2,Mb,T0],options);

% Save results
save x;
lengthx=size(x); Xinit=x(lengthx(1),:)';
save Xinit;
save C:/4s/AIChE//Ethanol_Benzene_Water/multivessel/decanting/Xinit.mat

Results=[Xinit(1+6+2*NT1+2*NT2), Xinit(2+6+3*NT1+2*NT2), Xinit(9+3*NT1+3*NT2); Xinit(1), Xinit(2+NT1), Xinit(3+NT1+NT2);...
         Xinit(1+3+NT1+NT2), Xinit(2+3+2*NT1+NT2), Xinit(6+2*NT1+2*NT2);...
         1-Xinit(1)-Xinit(3+NT1+NT2+1), 1-Xinit(2+NT1)-Xinit(3+NT1+NT2+2+NT1), 1-Xinit(3+NT1+NT2)-Xinit(6+2*NT1+2*NT2)];
