% programa 4.2, W.L.Luyben, pag. 98
% "example of iterative bubblepoint calculation
%  using Newton-Raphson algorithm"
% modif. E.Quintero.M 18nov'98
%------------------------------------------------------
% main program sets different values for initial
% guess of temperature
% specific chemical system is benzene/toluene 
% at 760 mm hg pressure
clear
format short e
% calculations of vapor pressure constants 
% for benzene and toluene
% data given at 100 and 125 degrees celcius
a1=log(2600/1360)/((1/(125+273))-(1/(100+273)));
b1=log(2600)-a1/(125+273);
a2=log(1140/550)/((1/(125+273))-(1/(100+273)));
b2=log(1140)-a2/(125+273);
% set liquid composition and pressure
x=0.5; p=760; t0=80; xx=[];
x,p
for nt=1:3
   t=t0;
   % subroutine bubpt
   loop=0;
   % temperature iteration loop
   while(1)
      loop=loop+1;
      if(loop>100) error('loop in bubpt subroutine'); end
      ps1=exp(a1/(t+273)+b1);
      ps2=exp(a2/(t+273)+b2);
      pcalc=x*ps1+(1-x)*ps2;
      % test for convergence
      if(abs(p-pcalc)>=p/10000)
         f=pcalc-p;
         df=-(x*ps1*a1+(1+x)*ps2*a2)/(t+273)^2;
         t=t-f/df;
      else
         y=x*ps1/p;
         break;
      end
   end
   xx=[xx;t0,t,y,loop];
   t0=t0+20;
end
xx