function [c,ceq]=f(w)
% This is file flashPHsrkeo.m 
%    Flash at given pressure and enthalpy
%    w = internal state vector of length 2*NC
% Usage for a case with NC=5 components (note: give composition only for NC-1 component):
%    x0= Initial estimates x, y0= Initial estimates yo; VF0= Initial estimate Vf; T0= Initial estimate T;
%    Zl0=initial estimate kiquid compressibility; Zg0= Initial estimate vapor compressibility; w0= [x0, y0, VF0, T0, Zl0, Zg0];    
%    Temperature=T, Vapor_fraction=VF, Liquid_Composition = x, Vapor_Composition = y
%Global variables for recursive routine
global hf p zf

% % Data (parameters and independent variables)
%Number of components
NC=5;
% zf % Specified feed composition 
% p % Specified pressure [N/m2]
% hf % Specified enthalpy [kJ/Kmol] (ref: h=0 for ideal gas at 298K)


% 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
T =w(2*NC);                            % Temperature [K]
Zl = w(2*NC+1);                        % Liquid Compressibility
Zg = w(2*NC+2);                        % Gas 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 = hf - VF*hg - (1-VF)*hl;                                  % = 0   Algebraic (1): Energy balance
f4 = [abs(errl) abs(errg)];  

c=[]; ceq=[f1';f2';f3';f4'];


