%Solve two-phase equilibria using Newton-Raphson iteration [x,y,l] = TpNewton(varargin). 
%The Hessian matrix is approximated  by ideal mixture theory. 
%The gradient canincorporate non-ideal behaviour.  Three cases are distinguished:
%
%             1) varargin = z,k
%             2) varargin = z,k,f,p
%             3) varargin = x,y,k,f,p,g,q
%
%z = feed composition vector
%k = K-value vector
%x = liquid composition vector (fractions on output)
%y = vapour composition vector (fractions on output)
%f = liquid activity or equation-of-state model
%p = liquid parameter struct
%g = vapor activity or equation-of-state model
%q = vapor parameter struct
%l = liquid fraction (mole liquid per mole feed)
%
%Copyright Tore Haug-Warberg 2003 (course SIK3035, NTNU)
%
 function [x,y,l] = TpNewton(varargin)

z  = varargin{1} + varargin{2};    % feed composition
nl = varargin{1}/sum(z); % normalized liq composition
nv = varargin{2}/sum(z); % normalized vap composition
z  = z./sum(z);         % normalized feed composition
k  = varargin{3};                          % K-values
f  = varargin{4};     % act.mod. or equation-of-state
p  = varargin{5};                  % model parameters
g  = varargin{6};     % act.mod. or equation-of-state
q  = varargin{7};                  % model parameters

tol = 1.e-6;                    % convergence tolerance
dm  = ones(length(k),1);        % chem. pot. 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; x=0.0*z; y=z; return; end
if max(k) < 1.0; l=1.0; y=0.0*z; x=z; return; end

while norm(dm)>tol & s>tol
x   = nl/sum(nl);
y   = nv/sum(nv);
dm  = log(y./x./k) + feval(g,y,q) - feval(f,x,p);
d   = (x.*y)./z;
dn  =-sum(nv)*sum(nl)*(d.*dm+d*(d'*dm)/(1-sum(d)));
s   = min(0.8/abs(min(min(dn./nv),min(-dn./nl))),1);
nv  = nv + s*dn;
nl  = z - nv;
end
l = sum(nl);

return