% running distillation model 'rectifier' for a ternary mixture
load Xinit                                   % loads state vector from the start_up period
% Number of stages in the rectifying section: 
Nr=51;
% Total number of 'stages'  
NT=Nr+2;

% Initial conditions
zFA0=Xinit(1:NT);                             % Initial composition of light component
zFB0=Xinit(NT+1:2*NT);                        % Initial composition of intermediate component 
M0=Xinit(2*NT+1:3*NT);                        % Initial molar holdup in each stage (kmol) 
T0=Xinit(3*NT+1:4*NT);                        % Initial temperature in each stage (K)

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


% Running column until steady state
[t,x]=ode15s(@input_rectifier,[0 3],[zFA0;zFB0;M0;T0],options);

% Save results
save x; 
lengthx=size(x); Xss=x(lengthx(1),:)';
save Xss;
Results=[Xss(2*NT+1), Xss(3*NT);...
         Xss(1), Xss(NT);...
         Xss(NT+1), Xss(2*NT);...
         1-Xss(1)-Xss(NT+1), 1-Xss(NT)-Xss(2*NT)];
