
function [c,ceq]=f(w)
% This is file flashTPsrkeo.m 
%    Flash at given pressure and temperature with explicit fugacities
%    w = internal state vector of length 2*NC-1
% Usage for a case with NC=5 components (note: give compositions only for NC-1 components):
%   x0=Initial estimate x, y0= initial estimate y; VF0= initial estimate vf; Zl0=initial estimate kiquid compressibility;
%   Zg0= Initial estimate vapor compressibility; w0= [x0, y0, VF0,Zl0,Zg0]; 
% Definition of global variables for recursive routine
global T p zf
% All in SI units 
% T = Temperature [K]
% p = Pressure [N/m^2]
% zf = Initial composition

%Number of components
NC=5;



% Extract states:
x(1:NC-1)=w(1:NC-1);                   % liquid composition 
y(1:NC-1)=w(NC:2*NC-2);                % vapor composition
VF =w(2*NC-1);                         % vapor fraction       
Zl =w(2*NC);                           % Liquid compressibility
Zg = w(2*NC+1);                        % Vapor compressibility

% Composition of last component:
x(NC)=1-sum(x(1:NC-1));        
y(NC)=1-sum(y(1:NC-1));                                              

% Fugacity, enthalpy and volumes from SRK
[phil,hl,Vlm,errl]=srkseo(x,T,p,Zl);       % liquid
[phig,hg,Vgm,errg]=srkseo(y,T,p,Zg);        % vapor

K = phil./phig;


f1 = y - K.*x;                                                % = 0   Algebraic (NC): VLE
f2(1:NC-1) = zf(1:NC-1) - VF.*y(1:NC-1) - (1-VF).*x(1:NC-1);  % = 0   Algebraic(NC-1): Component mass balances
f3 = [abs(errl) abs(errg)];  
c=[]; ceq=[f1';f2';f3'];

 