% running distillation model 'rectifier' for a ternary mixture

% Number of stages in the rectifying section: 
Nr = 50;
% Total number of 'stages'  
NT = Nr + 2;                                  % We include condenser and reboiler 
    
% Initial conditions
zFA=1/3*ones(1,NT);                           % Initial composition of light component
zFB=1/3*ones(1,NT);                           % Initial composition of intermediate component
Md=0.035;                                     % Initial holdup in condenser (kmol)
M1=1/500*ones(1,Nr);                          % Constant holdup in rectifying section (kmol)  
Mb=5.250;                                     % Initial holdup in reboiler (kmol) 
T0=350*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_rectifier,[0 20000],[zFA,zFB,Md,M1,Mb,T0],options);

% Save results
save x;
lengthx=size(x); 

% 'Xinit' is the matrix of the final state variables
% First NT states are molar fractions of the light component A, 
% Next NT states are molar fractions of the intermediate component B 
% (Molar fractions of the heavy component C are calculated from XC=1-XA-XB) 
% Next NT states are liquid holdups in the stages (including condenser + reboiler) 
% Next NT states are temperatures in the stages (including condenser + reboiler) 
Xinit=x(lengthx(1),:)';
save Xinit