function [c,ceq]=f(w)
% This is file flashPTsrk.m 
%    Flash at given pressure and temperature
%    w = internal state vector of length 2*NC-1
% Usage for a case with NC=4 components (note: give compositions only for NC-1 components):
%    x0=[0.3 0.4 0.1], y0=[0.52 0.36 0.07]; VF0=0.5; w0= [x0, y0, VF0];    % initial estimates of x, y and V/F
%    w=fmincon('1',w0,[],[],[],[],[],[],@flashPTsrk);
%    NC=(length(w)+1)/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);
%   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 
T=300;                                 % Specified temperature [K]    
p=5.e5;                                % Specified pressure [N/m2]

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              
% 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
  
c=[]; ceq=[f1';f2'];
 


