% This file calls the stiff differential equation routine ode15s, 
% to solve the distillation column problem in section 5.6;
%

 global coeff T xb xt xz JT right_hand;

tt=cputime;

% initial conditions
%
xb(1) = 0.3;	xt(1) = 0.4;
xb(2) = 0.25;	xt(2) = 0.35;
xb(3) = 0.25;	xt(3) = 0.35;
xb(4) = 0.25;	xt(4) = 0.35;
xb(5) = 0.2;	xt(5) = 0.3;


P = 1.0;	% atm
%
% parameters used to calculate HL and HV
%
MWb = 78.1; MWt = 92.1; MWz = 106.2;
lambda = 155;

% guess of initial T
%
for J = 1:5,
    T(J) = 200;
end

%
% matrix initialization
%
coeff = zeros(10);
right_hand = zeros(10, 1);

%
% parameters for matrix inversion
%
Q1 = 85000;
F = 2.5;
reflux = 2.0;
xtf = 0.30;
TF = 250;
xbf = (P - feval('Pz', TF) - xtf*(feval('Pt', TF)-feval('Pz', TF)))...
       /(feval('Pb', TF)-feval('Pz', TF));
HF = F*0.5*TF*(MWb*xbf + MWt*xtf + MWz*(1-xbf-xtf));

    for J = 1:5,
        xz(J) = 1 - xb(J) - xt(J);
     end

    %
    % matrix to solve for L and V
    %
    coeff(1,1) = 1.0;
    coeff(1,2) = -reflux;
    coeff(2,1) = -1.0;
    coeff(2,2) = -1.0;
    coeff(2,4) = 1.0;
    coeff(4,1) = 1.0;
    coeff(4,3) = -1.0;
    coeff(4,4) = -1.0;
    coeff(4,6) = 1.0;
    coeff(6,3) = 1.0;
    coeff(6,5) = -1.0;
    coeff(6,6) = -1.0;
    coeff(6,8) = 1.0;
    coeff(8,5) = 1.0;
    coeff(8,7) = -1.0;
    coeff(8,8) = -1.0;
    coeff(8,10) = 1.0;
    coeff(10,7) = 1.0;
    coeff(10,9) = -1.0;
    coeff(10,10) = -1.0;

    right_hand(5) = -HF;
    right_hand(6) = -F;
    right_hand(9) = -Q1;

    c0 = [xb(1) xt(1) xb(2) xt(2) xb(3) xt(3) xb(4) xt(4) xb(5) xt(5)]';

    [t, c] = ode15s('model56', [0, 0.8], c0);
%
% set up plotting vectors
%
    N = length(t);
for I = 1:N,
  
    xb(1)=c(I,1); xt(1)=c(I,2);
    xb(2)=c(I,3); xt(2)=c(I,4);
    xb(3)=c(I,5); xt(3)=c(I,6);
    xb(4)=c(I,7); xt(4)=c(I,8);
    xb(5)=c(I,9); xt(5)=c(I,10);

    xbplot(I) = xb(5);
    xtplot(I) = xt(5);
    xzplot(I) = 1-xb(5)-xt(5);

 for J = 1:5,
   JT = J;
   xz(J) = 1-xb(J)-xt(J);
   T(J) = fzero('impT',T(J));
   yb(J) = feval('Pb',T(J))/P*xb(J);
   yt(J) = feval('Pt',T(J))/P*xt(J);
   yz(J) = feval('Pz',T(J))/P*xz(J);
 end
ybplot(I) = yb(1);
ytplot(I) = yt(1);
yzplot(I) = yz(1);

  for J = 1:5,
    HL(J) = (MWb*xb(J)+MWt*xt(J)+MWz*xz(J))*(.5*T(J));
    HV(J) = (MWb*yb(J)+MWt*yt(J)+MWz*yz(J))*(.5*T(J)+lambda);
  end
coeff(3,1) = HL(1);
coeff(3,3) = -HL(2);
coeff(3,4) = -HV(2);
coeff(3,6) = HV(3);
coeff(5,3) = HL(2);
coeff(5,5) = -HL(3);
coeff(5,6) = -HV(3);
coeff(5,8) = HV(4);
coeff(7,5) = HL(3);
coeff(7,7) = -HL(4);
coeff(7,8) = -HV(4);
coeff(7,10) = HV(5);
coeff(9,7) = HL(4);
coeff(9,9) = -HL(5);
coeff(9,10) = -HV(5);
LV = coeff\right_hand;
L(1)=LV(1); L(2)=LV(3); L(3)=LV(5); L(4)=LV(7); L(5)=LV(9);
V(1)=LV(2); V(2)=LV(4); V(3)=LV(6); V(4)=LV(8); V(5)=LV(10);
v1plot(I) = V(1);
l5plot(I) = L(5);

end
time = cputime-tt

plot(t, v1plot);
title('Vapor Overhead Flow Rate (Example 5.6)');
xlabel('Time (hr)');
ylabel('Flow Rate (lb.mol/hr)');
%print -deps fig530.eps;
pause;
plot(t, l5plot);
title('Bottom Liquid Flow Rate')
xlabel('Time (hr)');
ylabel('Flow Rate (lb.mol/hr)');
%print -deps fig531.eps;
pause;
plot(t, ybplot, '-', t, ytplot, '--', t, yzplot, '*');
title('Overhead Vapor Composition');
xlabel('Time (hr)');
ylabel('Mole Fraction');
text(.04,.75,'benzene');
text(.32,.65,'toluene');
text(.4,.3,'xylene');
%print -deps fig532.eps;
pause;
plot(t, xbplot, '-', t, xtplot, '--', t, xzplot, '*');
title('Bottom Flow Composition');
xlabel('Time (hr)');
ylabel('Mole Fraction');
text(.07,.1,'benzene');
text(.2,.25,'toluene');
text(.3,.8,'xylene');
%print -deps fig533.eps;
