function [c,ceq]=f(w)
% This is file flashPHsrk.m 
%    Flash at given pressure and enthalpy
%    w = internal state vector of length 2*NC
% Usage for a case with NC=4 components (note: give composition only for NC-1 component):
%    x0=[0.5 0.5 0], y0=x0; VF0=0.5; T0=300; w0= [x0, y0, VF0, T0];    % initial estimates of x, y, V/F and T
%    w=fmincon('1',w0,[],[],[],[],[],[],@flashPHsrk);
%    NC=length(w)/2; x =  [w(1:NC-1) 1-sum(w(1:NC-1))]; y =  [w(NC:2*NC-2) 1-sum(w(NC:2*NC-2))]; VF=w(2*NC-1); T = w(2*NC);
%    Temperature=T, Vapor_fraction=VF, Liquid_Composition = x, Vapor_Composition = y

% All in SI units 
% Data (parameters and independent variables)
NC=4;
zf=[0.4 0.4 0.1 0.1];                  % Specified feed composition 
p=5.e5;                                % Specified pressure [N/m2]
hf = -10.e3;                           % Specified enthalpy [J/mol] (ref: h=0 for ideal gas at 298K)

liquid=1;vapor=2;
R=8.13;                                % J/mol K

% 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]

% 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
[Zl,phil,hl,Vlm]=srksis(x,T,p,liquid);       % liquid
[Zg,phig,hg,Vgm]=srksis(y,T,p,vapor);        % 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
  
c=[]; ceq=[f1';f2';f3'];


