function VLEfitting
% 估計二元混合物的van Laar模型中的參數
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/04/09 $
%
%   [Ref]: Edgar, T. F., Himmelblau D. M. and Lasdon L. S., Optimization 
%   of Chemical Processes, 2nd edition, McGraw-Hill, 2001（p.451）
%
%   x1,x2 = the liquid phase mole fraction of component 1 (water) and 2 (dioxane)
%   Pexp = the experimental data of total pressure, mmHg
%   T = temperature,℃
%   a1,a2 (vector) = Antoine constants of water and dioxane,respectively
%   Ps1, Ps2 = the saturation pressure of component 1 (water) and 2 (dioxane),mmHg

clear all
clc
x1 = 0.00:0.10:1.00; x2 = 1 - x1;
Pexp = [28.10,34.40,36.70,36.90,36.80,36.70,36.50,35.40,32.90,27.70,17.50];
T = 20;
a1 = [8.07131,1730.630,233.426];
a2 = [7.43155,1554.679,240.337];
Ps1 = Psat(T,a1);
Ps2 = Psat(T,a2);

% 非線性最小二乘法(非線性數據擬合)
lb = [0  0];
ub = [+inf  +inf];
beta0 = [1  1];
[beta,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,beta0,lb,ub,[],x1,x2,Ps1,Ps2,Pexp)           
ci = nlparci(beta,residual,jacobian)
[P1,P2,Pcal] = PressureCal(beta,x1,x2,Ps1,Ps2);
y1 = P1./Pcal;

% 結果
fprintf('Estimated Parameters:\n')
fprintf('\tA12 = %.3f ± %.3f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tA21 = %.3f ± %.3f\n',beta(2),ci(2,2)-beta(2))
fprintf('\tThe sum of the squares is: %.3f\n',resnorm)
fprintf('\tThe mean ΔP is: %.3f\n',mean(residual))
fprintf('The vapor phase mole fraction y1 is:\n')
fprintf('\t%.4f',y1)

% Plot
xx1 = linspace(x1(1),x1(end),200);
xx2 = 1 - xx1;
[P1,P2,Pcal] = PressureCal(beta,xx1,xx2,Ps1,Ps2);
yy1 = P1./Pcal;
plot(xx1,yy1)
xlabel('x_1')
ylabel('y_1')
refline([1,0])


% ------------------------------------------------------------------
function f = ObjFunc(A,x1,x2,Ps1,Ps2,Pexp)
[P1,P2,Pcal] = PressureCal(A,x1,x2,Ps1,Ps2);
f = Pcal - Pexp;

% ------------------------------------------------------------------
function [P1,P2,Pcal] = PressureCal(A,x1,x2,Ps1,Ps2)
gamma1 = exp(A(1)*(A(2)*x2./(A(1)*x1+A(2)*x2)).^2);
gamma2 = exp(A(2)*(A(1)*x1./(A(1)*x1+A(2)*x2)).^2);
P1 = x1.*gamma1*Ps1;
P2 = x2.*gamma2*Ps2;
Pcal = P1 + P2;

% ------------------------------------------------------------------
function Ps = Psat(T,a)     
% Calculation of the saturation pressure of component i
% T = temperature,℃
% a (vector) = Antoine constants of component i
logPs = a(1)-a(2)/(T+a(3));
Ps = exp(log(10)*logPs);