
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=5 components (note: give compositions only for NC-1 components):
%    x0=Initial estimate x, y0= initial estimate y; VF0= initial estimate vf;  w0= [x0, y0, VF0];    
%    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
% 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;

liquid=1;
vapor=2;
% 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]=srks(x,T,p,liquid);       % liquid
[Zg,phig,hg,Vgm]=srks(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'];

 