%This script calculate Homework 3 Pb 2
%created by Kyaw Tun
%Aug 15, 2001

%cgs unit used
%x1 = mole fration of liquid methane
%x2 = mole fration of liquid methanol
%y1 = mole fration of vapour methane
%x1 = mole fration of vapour methanol

global T;
global P;
global a;
global b;
global fid;

%first iteration assumtion
y1 = 1;
y2 = 0;

v1e = 36.0 ;			%partial molar volume of methane at infinate dilution
Pt = 20;					%total pressure
P = Pt;
P2s = 0.0401;			%vapour pressure of methanol
H12 = 1022;				%Hencry's constant
v2L = 0.792; 			%molor volume of methanol liquid
B11 = -53.9;
B12 = -166;
B22 = -4068;
T = 273.15;
R = 83.1451;

[fid,msg] = fopen('D:\D1.txt','w');
if fid == -1
   disp(msg);
   return;
end

phi1 = phiMe4;

x1 = y1 * phi1 * P / H12;
phi2s = phiMeOH2;

Bmix = y1^2 * B11 + 2*y1*y2 * B12 + y2^2 * B22;
phi2 = exp((2*y1*B12 - Bmix) * Pt/R/T);
y2 = (1-x1) * P2s * phi2s * exp(v2L * (Pt - P2s) /R/T) / phi2/Pt;
y1 = 1 - y2;

phi2 = .9;

fprintf(fid,'Homework 3: Set D No.1 \n\nRecord of iterations\n\n');
fprintf(fid,'y2\t\t   phi1\t\t   x1\t\t\t   phi2\n');
fprintf(fid,'0\t\t   %1.9f\t   %1.9f\t   %1.9f\n', phi1, x1, phi2);

for i=1:30
   Bmix = y1^2 * B11 + 2*y1*y2 * B12 + y2^2 * B22;
   phi1 = exp((2*y2*B12 - Bmix) * Pt/R/T);
   x1 = y1 * phi1 * Pt / (H12 * exp(v1e * (Pt - P2s) /R/T));
   phi2 = exp((2*y1*B12 - Bmix) * Pt/R/T);
   phi2s = exp((2*y1*B12 - Bmix) * P2s/R/T);
   y2 = (1-x1) * P2s * phi2s * exp(v2L * (Pt - P2s) /R/T) / phi2/Pt;
   y1 = 1 - y2;
   fprintf(fid,'%1.9f\t   %1.9f\t   %1.9f\t   %1.9f\n',y2, phi1, x1, phi2);
end

status = fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function phi = phiMe(v)
%Problem Set D1
%fugicity coefficient of methanol


format long;
global T;
global P;
global a;
global b;
global fid;

%SI unit used throught
Tc = 190.4;
Pc = 46 * 10^5;
omega = 0.011;
R = 8.31451;

%calculate Peng-robison constants
b = 0.07796 * R * Tc / Pc;
beta = 0.37464 + 1.54226 * omega - 0.26992 * omega^2;
alpha = (1 + beta * (1 - sqrt(T / Tc)))^2;
aTc = 0.457235 * (R * Tc)^2 / Pc;
a = aTc * alpha;
fprintf(fid,'\nProperties of Methane\na = %1.9g\nb = %1.9g\n', a, b);

   guess = R * T / P;
   v = fzero('PRvol', guess);
   z = P * v / (R * T);
   phi = exp(((a / (sqrt(8) * b)) * log(abs(((v+b)-sqrt(2)*b) / ((v+b)+sqrt(2)*b))) ...
      - R*T * log((v-b) / v) - R*T * (z - 1)) /(R*T));
 fprintf(fid,'v = %1.9g\nz = %1.9g\nfugacity coefficient of methane for the first iteration = %1.9g\n', v, z, phi);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%555
function phi = phiMe(v)
%Problem Set D1
%fugicity coefficient of methanol
%created by Kyaw Tun
%Aug 15, 2001

format long;
global T;
global P;
global a;
global b;
global fid;

%SI unit used throught
Tc = 190.4;
Pc = 46 * 10^5;
omega = 0.011;
R = 8.31451;

%calculate Peng-robison constants
b = 0.07796 * R * Tc / Pc;
beta = 0.37464 + 1.54226 * omega - 0.26992 * omega^2;
alpha = (1 + beta * (1 - sqrt(T / Tc)))^2;
aTc = 0.457235 * (R * Tc)^2 / Pc;
a = aTc * alpha;
fprintf(fid,'\nProperties of Methane\na = %1.9g\nb = %1.9g\n', a, b);

   guess = R * T / P;
   v = fzero('PRvol', guess);
   z = P * v / (R * T);
   phi = exp(((a / (sqrt(8) * b)) * log(abs(((v+b)-sqrt(2)*b) / ((v+b)+sqrt(2)*b))) ...
      - R*T * log((v-b) / v) - R*T * (z - 1)) /(R*T));
 fprintf(fid,'v = %1.9g\nz = %1.9g\nfugacity coefficient of methane for the first iteration = %1.9g\n', v, z, phi);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
function y = PRvol(v)
%This funciton calculate volume of gas using
%Peng-Robison equation of state


global P;
global T;
global a;
global b;
R = 8.3145;				%gas constant
y = P - R * T / (v - b) + a / (v * (v + b) + b * (v - b));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5