txy_diagram (AntCon,vLconst,P,isPlot)

function [xA,yA] = txy_diagram (AntCon,vLconst,P,isPlot)
%Calculate isobaric V-L equilibrium data for binary system
% for real mixture and van Laar equations
%
%   [xA,yA] = txy_diagram (AntCon,vLconst,P,isPlot)
%       AntCon  - antoine constants
%       vLconst - van Laar constants, form [A B]
%       P       - pressure
%       isPlot  - plot t-xy diagram (true/false)
%
%       xA - molar fraction of component A
%       yA - equilibrium molar fraction in Vapor phase

A = vLconst(1);
B = vLconst(2);

%mole fractions
xA = 0:0.01:1;
xB = 1-xA;

%activity coefficients
gA = [0, 10.^(A./((1+xA(2:end-1).*A./(1-xA(2:end-1))./B).^2)), 1];
gB = [1, 10.^(B./((1+(1-xA(2:end-1)).*B./xA(2:end-1)./A).^2)), 0];

%force column vector
xA = xA(:);
xB = xB(:);
gA = gA(:);
gB = gB(:);

%boiling point of pure compound
tA = AntCon(1,2)/(AntCon(1,1)-log10(P))-AntCon(1,3);
tB = AntCon(2,2)/(AntCon(2,1)-log10(P))-AntCon(2,3);

x0   = linspace(tB,tA,length(xA));
spec = optimset('Display','notify');
Temp = fsolve(@ucelovaFun,x0,spec,AntCon,xA,xB,gA,gB,P);
yA   = compute_yA(Temp,AntCon,xA,gA,P);

if isPlot
    h = figure;
    set(h,'Color',[1 1 1]);
    plot(xA,Temp,'b',yA,Temp,'r')
    xlabel('xA, yA')
    ylabel('T')
    title('T-x,y diagram')
    xlim([0 1])
end


function F = ucelovaFun (t,AC,xA,xB,gA,gB,P)
%nonlinear system of equations to solve
pA0 = antoineEq(AC(1,:),t);
pB0 = antoineEq(AC(2,:),t);
pA0 = pA0(:);
pB0 = pB0(:);

pA = pA0.*xA.*gA;
pB = pB0.*xB.*gB;
F  = pA+pB-P;

function F = compute_yA (t,AC,xA,gA,P)
% computes yA for xA
pA0 = antoineEq(AC(1,:),t);
pA0 = pA0(:);
F   = pA0/P.*xA.*gA;

function F = antoineEq (AConstants,t)
%Antoine equation
t = t(:);
F = 10.^(AConstants(1)-(AConstants(2)./(AConstants(3)+t)));