function [Xaq,Xorg,Laq] = LLE(z,temp)
% ********************************************************************************************************* %
% Copyright Tore Haug-Warberg 2003 (course SIK3035, NTNU)
% Solves liquid-liquid equilibria using Newton-Raphson iteration. 
% The Hessian matrix is approximated  by ideal mixture theory. 
% The gradient can incorporate non-ideal behaviour
% ******************************************************************************************************** %
% Program modified by Stathis Skouras, April 2003 
% ******************************************************************************************************** %
% System:                 (1) Ethanol  
%                         (2) Benzene  (organic phase)
%                         (3) Water    (aqueous phase)
% ********************************************************************************************************* %
% Gas constant
R=1.98721;   % universal gas constant [cal/molK]
% *************************************************************************
% Data that have to be changed for another mixture
Xaq=[0.35;0.05;0.6];            % An estimation of the distribution between the phases (stable up to the critical point)  
Xorg=[0.2;0.7;0.1];           % An estimation of the distribution between the phases (stable up to the critical point)
k  = ones(size(z));              % same reference state in both phases

% The NRTL activity model is used for both liquid phases
liq     = @ActivityNRTL;                                % function name

s.alpha = [0      0.2    0.2;                   % non-random par
           0.2    0      0.2;                   % from LLE data
           0.2    0.2    0];              

s.tau   = [  0.0      87.744  -441.74; ...         % interaction coefficient parameters [K]
           118.04      0.0       986.99; ...         % from LLE data
           376.33   2797.7         0.0];            

s.t     = 25+273.15;                                    % temperature [K]
% End of data that have to be changed
% ************************************************************************ %

% Parameters for the liquid-liquid algorithm
naq  = Xaq.*z;                  % Moles in aqueous phase
norg  = Xorg.*z;                % Moles in organic phase
f  = liq;                       % activity coefficient model for both phases (NRTL is used)
p  = s;                         % parameters for NRTL
g  = liq;                       % activity coefficient model for both phases (NRTL is used)
q  = s;                         % parameters for NRTL

tol = 1.e-6;                   % convergence tolerance
dm  = ones(length(k),1);       % chemical potensials difference
s   = 1;                       % relative step length

if min(z) < 0.0;                      return; end
if min(k) < 0.0;                      return; end
if min(k) > 1.0; l=0.0; Xaq=0.0*z; Xorg=z; return; end
if max(k) < 1.0; l=1.0; Xorg=0.0*z; Xaq=z; return; end

while norm(dm)>tol & s>tol
  Xaq   = naq/sum(naq);
  Xorg   = norg/sum(norg);
  dm  = log(Xorg./Xaq./k) + feval(g,Xorg,q) - feval(f,Xaq,p);
  d   = (Xaq.*Xorg)./z;
  dn  =-sum(norg)*sum(naq)*(d.*dm+d*(d'*dm)/(1-sum(d)));
  s   = min(0.8/abs(min(min(dn./norg),min(-dn./naq))),1);
  norg  = norg + s*dn;
  naq  = z - norg;
end
Laq = sum(naq);
return;