function [y,x,N]=McCabeT(alpha,R,xd,xb,zf,q)
% This function uses analytical solution for binary system
% alpha - relative volatility
% R     - reflux ratio
% xd    - distillate comp. mole fraction
% xb    - bottom comp. mole fraction
% zf    - feed composition
% q     - feed thermal condition

%E.g. [y,x,N]=McCabeT(3.8,1.2,0.95,0.05,0.45,0.55)

%Solve for (xi,yi) and store in vector xi
%i.e. intersection of q-line with enriching section 

if (q==1)
    xi(1)=zf; xi(2)=(R/(R+1))*zf+xd/(R+1);
else
    a(1,1) = R/(R+1); a(1,2) = -1; b(1) = -xd/(R+1);
    a(2,1) = q/(q-1); a(2,2) = -1; b(2) = zf/(q-1);
    xi=a\b'
end

x = xd; xx=x; xold=xd;
yy = xd; %Total condenser
y=yy;

D_F = (zf-xb)/(xd-xb);

L_V   = R/(R+1);                       %enriching section slope

% L_V_s = (R*D_F + q)/((R+1)*D_F-(1-q)) Stripping section slope Dr Kumar

L_V_s = (xi(2)-xb)/(xi(1)-xb);          %Stripping section slope
E_i=xd/(R+1);                          %intercept of enriching section
S_i=(1-L_V_s)*xb;                      %intercept of stripping section

xx = Ialpha(yy,alpha); x = [x,xx];     % equili curve

yp(1)=xd;xp(1)=xd;                     % first two points, the second one
xp(2)=xx;yp(2)=yp(1);                  % on the equilibrium curve

while (xx >= xb)                       % keep doing calc all the way till
                                       % it reaches a composition smaller
                                       % than Xb
                                   
   if (xx >= xi(1)),                   % enriching section: is the X 
                                       % estimated from the alpha model
                                       % higher that the interception
                                       % composition?
                                       
                                       % ***** YES *****
                                       
     yy=L_V*xx + E_i; y=[y,yy];        % Passing stream composition op line
     
     yp=[yp,yy];xp=[xp,xx];            % store the points
     
     xx = Ialpha(yy,alpha); x = [x,xx];%equili curve
     
     yp=[yp,yy];xp=[xp,xx];            % store the points
     
     if (abs(xx-xold) <= 1.e-7)
        error(' Possible pinch')       % pinch???
     end
     
     xold = xx;                        % store this as the last one calculated
     
                                       % ***** NO *****
     
  else                                 % the next point
                                       % will rest on the stripping section
     yy=L_V_s*xx + S_i; y=[y,yy];      % op line
          
     yp=[yp,yy];xp=[xp,xx];            % store the points
     
     xx = Ialpha(yy,alpha); x = [x,xx];%equilibrium 
     yp=[yp,yy];xp=[xp,xx];
  end
end

N=length(x)-1;                         %count number of stages, it removes
                                       % one because it's the first point
                                       % (xd, xd) is not a step
clf;

%Graphics construction begins
axis([0 1 0 1])
%45 deg line
x45(1)=0;y45(1)=0;x45(2)=1;y45(2)=1;
plot(x45,y45,'k'); hold on

%xd
xt(1)=xd, yt(1)=0; xt(2)=xd;yt(2)=xd; plot(xt,yt,'r');

%xb
xt(1)=xb, yt(1)=0; xt(2)=xb;yt(2)=xb; plot(xt,yt,'r');

%zf
xt(1)=zf, yt(1)=0; xt(2)=zf;yt(2)=zf; plot(xt,yt,'r');

%q-line
xt(1)=zf;yt(1)=zf; xt(2)=xi(1);yt(2)=xi(2);plot(xt,yt);

% Enriching section operating line
xt(1)=xd; yt(1)=xd;
xt(2)=0; yt(2)=E_i;
plot(xt,yt)

%Equilibrium curve
xeq=0:.05:1;
yeq=Eq_alpha(xeq,alpha);
plot(xeq,yeq);

%Staircase construction: all the points stored
plot(xp,yp)

% Stripping section operating line
xs(1)=xi(1);ys(1)=xi(2);
xs(2)=xb;ys(2)=xb;
plot(xs,ys);

xlabel('x'); ylabel('y');

ans=input('Want a composition plot?','s')
if (ans == 'y')
    clf;
    plot(1:N,x(1:N),'o',1:N,y(1:N),'+');
    xlabel('Stage number'); ylabel('Composition');
end
%Graphics construction ends

function y=Eq_alpha(x,alpha)
 y=alpha*x./(1+(alpha-1)*x);
return

function x=Ialpha(y,alpha)
 x=y/(alpha-(alpha-1)*y);
return

