function [c,ceq]=f(w)
% This is file flashPTsrk1.m 
% -----------------------------------------------------------------
% MODIFIED MARCH 2014 BY S. SKOGESTAD: Includes zl and zv AS STATES (thus avoid solving the cubic SRK-equation to find z)
% The modification is very efficient for dynami simulation, but it requires good initial estimates + SCALING needs to be reasonable (which it is not in this version)
% -----------------------------------------------------------------
%    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):
    NC=4; x0=[0.31 0.44 0.12], y0=[0.57 0.33 0.07]; VF0=0.35; Zl0=0.02; Zg0=0.9; w0= [x0, y0, VF0, Zl0, Zg0];    % initial estimates of x, y and V/F + zl, zv
    UB(1:2*NC+1)=1; UB(2*NC+1)=3;
    LB(1:2*NC+1)=0;
    w=fmincon('1',w0,[],[],[],[],LB,UB,@flashPTsrk1);
    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]

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
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  + error in cubic EOS from SRK

[phil,hl,Vlm,errl]=srksis1(x,T,p,Zl);      % liquid
[phig,hg,Vgm,errg]=srksis1(y,T,p,Zg);       % vapor

errl
errg

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)];                                            % = 0   Algebraic (2): Cubic equation for each phase
  
c=[]; ceq=[f1';f2';f3'];
 


