function [xs,ys,vof]=flashr(t,p,zs,erc)
% flashr                - Performs a flash calculation with a table of results
% function [xs,ys,vof]=flashr(t,p,zs,erc)
% Argument List
% Argument   Gives
%    t      temperature in Tdeg.
%    p      pressure in kPa.
%    zs     vector of feed mol fractions.
%    erc    gives fzero the difference between successive values of x
%           allowed when convergence is deemed to be achieved.
%  If erc is not given, then the program uses erc = .0001
%
%    xs is the returned vector of liquid phase component fractions
%    ys is the returned vector of gas phase component fractions
%    vof is the mass fraction gas in original feed
% Example:>> start301
%         >> [xs ys vof] = flashr(375,101,[.4 .4 .2]);
%         >> z = [xs;ys];
%         >> bar(z);
%  This will produce a bar graph displaying the relative
%  compositions of the liquid (xs) and the vapor (ys) phases
global cnms Tdeg
if nargin<=3
   erc=.0001;
end
ks=kval(t,p);
omks=1-ks;
P1=[zs;omks];
fx='sum(P1(1,:).*P1(2,:)./(1-x.*P1(2,:)))';
x=0;
f1=eval(fx);
x=1;
f2=eval(fx);
if f1*f2>=0
   disp('There is no V/F for this flash.')
   return
end
vof=fzero(fx,[0 1],erc,[],P1);
xs=zs./(1-vof.*c(2,:));
fprintf(['Temperature= %6.2f',Tdeg,' Pressure= %6.2f kPa V/F=%6.4f\n'],t,p,vof)
%disp('Temperature  Pressure    V/F')
%disp([t,p,vof])
disp('Compound   Feed      Liquid    Vapor     K=y/x')
ys=xs.*ks;
m=[zs;xs;ys;ks]';
cm=charmx(m,10,5);
for k=1:length(zs)
   disp([cnms(k,:),cm(k,:)])
end
