%octave
1;
function psat=psat(a,b,c,t)
    lnpsat=a-b/(t+c);
    psat=exp(lnpsat);
endfunction
    
function tsat=tsat(a,b,c,p)
    tsat=(b/(a-log(p)))-c;
endfunction

% Presión de burbuja
% Inputs:
x=[0.05 0.95];
T=343.15-273.15; % ºC

% Coeficientes de actividad:
gamma1=exp(0.95*x(2)^2);
gamma2=exp(0.95*x(1)^2);
gammas=[gamma1 gamma2];
% Cálculo de la presión:
psats=[79.8 40.50];
P=sum(x.*psats.*gammas)
y=(x.*psats.*gammas)/P


%Presión de rocío:
y=[0.05 0.95];
T=343.15;
psats=[79.8 40.5];
gammas=[1 1];
Pguess=1/sum(y./(psats.*gammas));
tol=0.001;
while 1
    x=(y*Pguess)./(psats.*gammas);
    gammasn=[exp(0.95*x(2)^2) exp(0.95*x(1)^2)];
    if abs(gammasn-gammas)<=tol
        break
    else
        gammas=gammasn;
    end
end
P=1/sum(y./(psats.*gammas))
x=(y*P)./(psats.*gammas)

%Temperatura de burbuja:
x=[0.5 0.5];
P=101.325;
tsats=[tsat(14.3145,2756.22,228.060,P) tsat(16.5785,3638.27,239.5,P)];
tguess=mean(tsats);
tol=0.001;
while 1
    gammas=[exp(0.605*x(2)^2) exp(0.605*x(1)^2)];
    psats=[psat(14.3145,2756.22,228.060,tguess) psat(16.5785,3638.27,239.5,tguess)];
    alfa=psats(1)/psats(2);
    psat2=P/(x(1)*alfa*gammas(1)+x(2)*gammas(2));
    tsat2=tsat(16.5785,3638.27,239.5,psat2);
    if abs(tguess-tsat2)<=tol
        break
    else
        tguess=tsat2;
    end
end
T=tsat2
y=(x.*psats.*gammas)/P

%Temperatura de rocío:
y=[0.5 0.5];
P=101.325;
gammag=[1 1];
tsats=[tsat(14.3145,2756.22,228.060,P) tsat(16.5785,3638.27,239.5,P)];
tguess=mean(tsats);
tol=0.001;
while 1
    psats=[psat(14.3145,2756.22,228.060,tguess) psat(16.5785,3638.27,239.5,tguess)];
    alfa=psats(1)/psats(2);
    x=(y*P)./(gammag.*psats);
    gammasn=[exp(0.605*x(2)^2) exp(0.605*x(1)^2)];
    psat1=P*((y(1)/gammasn(1))+(y(2)/(gammasn(2)))*alfa);
    tsat1=tsat(14.3145,2756.22,228.060,psat1);
    if abs(tguess-tsat1)<=tol
        break
    else
        tguess=tsat1;
    end
end
T=tsat1
x=(y*P)./(gammasn.*psats)

% Presión de burbuja
% Inputs:
x=[0.65 0.35];
%T=; % ºC

% Coeficientes de actividad:
gamma1=exp(0.67*x(2)^2);
gamma2=exp(0.67*x(1)^2);
gammas=[gamma1 gamma2];
% Cálculo de la presión:
psats=[32.27 73.14];
P=sum(x.*psats.*gammas)
y=(x.*psats.*gammas)/P


%Presión de rocío:
y=[0.65 0.35];
%T=343.15;
psats=[32.27 73.14];
gammas=[1 1];
Pguess=1/sum(y./(psats.*gammas));
tol=0.001;
while 1
    x=(y*Pguess)./(psats.*gammas);
    gammasn=[exp(0.95*x(2)^2) exp(0.95*x(1)^2)];
    if abs(gammasn-gammas)<=tol
        break
    else
        gammas=gammasn;
    end
end
P=1/sum(y./(psats.*gammas))
x=(y*P)./(psats.*gammas)

z=[0.65 0.35];
k=[0.6485 2.0547];
vlg=0.5;
tol=0.001;
while 1
    lx=z./(1+vlg*k);
    l=sum(lx);
    v=1-l;
    vln=v/l;
    if abs(vlg-vln)<=tol
        break
    else
        vlg=vln;
    end
end
v
x=lx/l