function [x, y, VF, T  ] = flashPH( P, h,  x0, y0, z, T0, VF0 )
% %Recursive function for the claculation of a PH flash.
% Extracts liquid composition, vapor composition, vapour
% fraction and temperature for liquid and for vapor phase, given Pressure and enthalpy. 
% Definition of global variables for recursive routine.
global p zf hf
p=P; hf=h; zf=z;
% zf % Specified feed composition 
% p % Specified pressure [N/m2]
% hf % Specified enthalpy [kJ/Kmol] (ref: h=0 for ideal gas at 298K)

% Initialization vector 
w0= [x0, y0, VF0, T0]; 
%x0 initial estimate for liquid composition
%y0 initial estimate for vapor composition
%vf0 initial estimate for vaporr fraction
%T0 intial estimate for temperature

%Boundary constrains. All the states except temperature must lie between 0 and 1. 
lb=zeros(size(w0));
ub=ones(size(w0));
ub(10)=Inf;
%Inequality constrains. The added compositions cannot be larger than 1.
A=[1 1 1 1 0 0 0 0 0 0;
    0 0 0 0 1 1 1 1 0 0];
B=[1; 1];
%Solver options
options = optimoptions('fmincon','Algorithm','interior-point','TolCon', 1e-6,'TolFun', 1e-4, 'MaxIter', 3000, 'TolX', 1e-12, 'ScaleProblem','obj-and-constr');

%Solution subject to the equilibrium and solving the EOS. For further
%details check flashPHsrk.m
[w,fval,exitflag,output]=fmincon('1',w0,A,B,[],[],lb,ub,@flashPHsrk,options);

% %Data extraction
NC=(length(w)+1)/2; 

%Liquid and vapor compositions
x =  [w(1:4) 1-sum(w(1:4))];
y =[w(5:8) 1-sum(w(5:8))];
%Vapor fraction
VF=w(9); 
%Temperature
T = w(10);
end

