%Experimental  data  for  Hexane(1) - Toluen(2)  taken  from  the 
%"Dechema Data Series". Temperatures are in Celsius, pressures in
%mmHg.
 t  = [105.4 104.4 103.1 101.7 100.3  98.0  97.8 95.8...
        93.6  92.4  91.6  89.2  85.0  82.8  80.0 78.3...
        77.0  76.3  74.9  72.7  71.4  70.2];
 x1 = [.047  .057  .071  .087  .104  .125  .141 .172 ...
       .201  .224  .240  .290  .379  .439  .516 .582 ...
       .639  .659  .754  .822  .885  .938];
 y1 = [.154  .176  .212  .251  .301  .331  .364 .424 ...
       .471  .512  .531  .592  .667  .717  .766 .812 ...
       .833  .851  .886  .921  .954  .975];
%
 p  = 760;                                      % system pressure
 t  = [105.4 101.7 97.8 92.4 85.0 78.3 74.9 70.2];      % Celsius
 x1 = [.047  .087  .141 .224 .379 .582 .754 .938];       % liquid
 y1 = [.154  .251  .364 .512 .667 .812 .886 .975];        % vapor
% 
 c1 = [6.91058 1189.64 226.280];       % published Antoine Hexane
 c2 = [6.95087 1342.31 219.187];      % published Antoine Toluene
 t1 = c1(2)/(c1(1)-log10(p))-c1(3);                  % nbp Hexane
 t2 = c2(2)/(c2(1)-log10(p))-c2(3);                 % nbp Toluene
 k  = [c1(1)-c1(2)./(t+c1(3));...          % log10(p/mmHg) Hexane
       c2(1)-c2(2)./(t+c2(3))];           % log10(p/mmHg) Toluene
 k  = 10.^k/p;                        % K-values Hexane & Toluene
 x  = [x1;1-x1];               % initialize liquid mole fractions
 y  = [y1;1-y1];               % initialize vapour mole fractions
% 
%Assume ideal liquid mixture and ideal gas phase.This is equival-
%ent to assuming Raoult's law (constant K-values)
 for i=1:length(t)
   [x(:,i),y(:,i)] = tpkvalue(x(:,i)+y(:,i),k(:,i));
 end
 ss = norm(x(1,:)-x1) + norm(y(1,:)-y1);
% 
%Show experimental points
 subplot('position',[.1 .12 .36 .80])
 plot(x1,t,'ob',y1,t,'vr'), hold on;
 axis([0,1,60,120])
%
%Show calculated results
 plot([0,x(1,:),1],[t2,t,t1],'--k');
 plot([0,y(1,:),1],[t2,t,t1],'--k');
 text(0.4,112,['ss1=',num2str(ss,3)],'BackgroundColor','w')
% 
%The van Laar activity model is used in the liquid phase
 liq   = @ActivityBinVanLaar;             % activity coefficients
 vap   = @ActivityRaoult;                         % ideal mixture
 s.a12 = 0.2222;                  % published model parameter A12
 s.a21 = 0.4306;                  % published model parameter A21
 for i=1:length(t)
%  [x(:,i),y(:,i)] = tpnewton(x(:,i),y(:,i),k(:,i),liq,s,vap,[]);
   [x(:,i),y(:,i)] = tpkvalue(x(:,i),y(:,i),k(:,i),liq,s,vap,[]);
 end
 ss = norm(x(1,:)-x1) + norm(y(1,:)-y1);
%
%Show calculated results
 plot([0,x(1,:),1],[t2,t,t1],'b');
 plot([0,y(1,:),1],[t2,t,t1],'r'), hold off;
 text(0.4,107,['ss2=',num2str(ss,3)],'BackgroundColor','w')
 xlabel('hexane')
 ylabel('Celsius')
 title('published')
% 
%The Antoine pressure is quite off for toluene. Redo the calcula-
%tions after adjusting  the  pressure 2.0 percent
 off    = 0.020;               % offset in Antoine vapor pressure
 c2(1)  = c2(1) + off;                 % adjusted Antoine Toluene
 t2     = c2(2)/(c2(1)-log10(p))-c2(3);             % nbp Toluene
 k(2,:) = k(2,:)*10^off;               % adjusted K-value Toluene
% 
%Raoults law.
 for i=1:length(t)
%  [x(:,i),y(:,i)] = tpnewton(x(:,i)+y(:,i),k(:,i));
   [x(:,i),y(:,i)] = tpkvalue(x(:,i)+y(:,i),k(:,i));
 end
 ss = norm(x(1,:)-x1) + norm(y(1,:)-y1);
% 
%Show experimental points
 subplot('position',[.62 .12 .36 .80])
 plot(x1,t,'ob',y1,t,'vr'), hold on;
 axis([0,1,60,120])
%
%Show calculated results
 plot([0,x(1,:),1],[t2,t,t1],'--k');
 plot([0,y(1,:),1],[t2,t,t1],'--k');
 text(0.4,112,['ss3=',num2str(ss,3)],'BackgroundColor','w')
%
%van Laar activity.
 for i=1:length(t)
%  [x(:,i),y(:,i)] = tpnewton(x(:,i),y(:,i),k(:,i),liq,s,vap,[]);
   [x(:,i),y(:,i)] = tpkvalue(x(:,i),y(:,i),k(:,i),liq,s,vap,[]);
 end
 ss = norm(x(1,:)-x1) + norm(y(1,:)-y1);
%
%Show calculated results
 plot([0,x(1,:),1],[t2,t,t1],'b');
 plot([0,y(1,:),1],[t2,t,t1],'r'), hold off;
 text(0.4,107,['ss4=',num2str(ss,3)],'BackgroundColor','w')
 xlabel('hexane')
 xlabel('hexane')
 title('corrected')
 tabbe('HexaneToluene',10,6,'Times',11,0.8,6,'on',':');
