function [c,ceq]=f(w)
% This is file bubsrk.m 
%    Bubble point flash at given temperature (can easily change to given pressure or dew point)
%    w = internal state vector of length NC
% Usage for a case with NC=4 components (note: give y0 only for NC-1 component):
%    y0=[0.5 0.5 0], p0=10.e5; w0= [y0, p0];     % initial estimates of y and p
%    w=fmincon('1',w0,[],[],[],[],[],[],@bubsrk);
%    NC=length(w); p=w(NC)/1.e5; y =  [w(1:NC-1) 1-sum(w(1:NC-1))]
%    Bubble_pressure_in_bar=p, Vapor_Composition = y

% All in SI units 
% Data (parameters and independent variables)
NC=4;
x=[0.5 0.5 0 0];                       % Specified liquid composition (if DEW: change x to y)
T=300;                                 % Specified temperature [K]    (if given p: change T to p)

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

% Extract states:
y(1:NC-1)=w(1:NC-1);                   % vapor composition            (if DEW: change y to x)
p =w(NC);                              % pressure [N/m2]              (if given p: change p to T)
% Composition of last component:
y(NC)=1-sum(y(1:NC-1));                %                              (if DEW: change y to x)

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

F = y - K.*x;                          % = 0	   Algebraic (NC): VLE 
  
c=[]; ceq=F';


