main_pr.m

clear all; clc

% Input (C3, n-C4, n-C5, n-C6)
p  = 5; % bar
T  = 95; % C
x  = [0.027361136 0.112926395 0.290338156 0.569374312]';
y  = [0.155276407 0.266260893 0.307352428 0.271110272]';
nc = length(x);

% Parameters
pc     = [616.410385 550.70829 487.326798 436.56359]'; % psia
pc     = pc*1e5/14.5038; % Pa
Tc     = [206.006002 305.618002 385.412001 453.47]'; % F
Tc     = (Tc - 32)/1.8 + 273.15; % K
w      = [0.152 0.201 0.254 0.301]';
k(1,1) =  0;
k(1,2) =  3.3000e-3;
k(1,3) =  0.0267;
k(1,4) =  7.0000e-4;
k(2,2) =  0;
k(2,3) =  0.0174;
k(2,4) = -5.6000e-3;
k(3,3) =  0;
k(3,4) =  0;
k(4,4) =  0;
k      =  k + k';
cpig   = [51.92 71.34 88.05 104.40; 192.45 243.00 301.10 352.30; 1626.50 1630.00 1650.20 1694.60;
    116.80 150.33 189.20 236.90; 723.60 730.42 747.60 761.60]'; % J/mol-K
DHf    = [-1.0468e+5 -1.2579e+5 -1.4676e+5 -1.6694e+5]'; % J/mol

tic
for i = 1:100
    [VL,ZL,phiL1,HL1] = pr(x,p*1e5,T+273.15,pc,Tc,w,k,cpig,DHf,'L');
    [VV,ZV,phiV1,HV1] = pr(y,p*1e5,T+273.15,pc,Tc,w,k,cpig,DHf,'V');
    K1 = phiL1./phiV1;
end
toc