VanDerWaalDensityCalculator

VDWDensityCalc.m

clear;
clc;


T = linspace(0,500,2000);
V1 = [];
V2 = [];
V3 = [];

R = 8.31446261815324E-2; 	#L.bar.K^-1.mol^-1
a = 1.382;
b = 0.03186;
P = 120.0; #Bar
M = 32.2;

for ii = 1:length(T)
  V1(ii) = -(-3*a/P + (-P*b - R*T(ii))^2/P^2)/(3*(sqrt(-4*(-3*a/P + (-P*b - R*T(ii))^2/P^2)^3 + (-27*a*b/P - 9*a*(-P*b - R*T(ii))/P^2 + 2*(-P*b - R*T(ii))^3/P^3)^2)/2 - 27*a*b/(2*P) - 9*a*(-P*b - R*T(ii))/(2*P^2) + (-P*b - R*T(ii))^3/P^3)^(1/3)) - (sqrt(-4*(-3*a/P + (-P*b - R*T(ii))^2/P^2)^3 + (-27*a*b/P - 9*a*(-P*b - R*T(ii))/P^2 + 2*(-P*b - R*T(ii))^3/P^3)^2)/2 - 27*a*b/(2*P) - 9*a*(-P*b - R*T(ii))/(2*P^2) + (-P*b - R*T(ii))^3/P^3)^(1/3)/3 - (-P*b - R*T(ii))/(3*P);
  V2(ii) = -(-3*a/P + (-P*b - R*T(ii))^2/P^2)/(3*(-1/2 - sqrt(3)*I/2)*(sqrt(-4*(-3*a/P + (-P*b - R*T(ii))^2/P^2)^3 + (-27*a*b/P - 9*a*(-P*b - R*T(ii))/P^2 + 2*(-P*b - R*T(ii))^3/P^3)^2)/2 - 27*a*b/(2*P) - 9*a*(-P*b - R*T(ii))/(2*P^2) + (-P*b - R*T(ii))^3/P^3)^(1/3)) - (-1/2 - sqrt(3)*I/2)*(sqrt(-4*(-3*a/P + (-P*b - R*T(ii))^2/P^2)^3 + (-27*a*b/P - 9*a*(-P*b - R*T(ii))/P^2 + 2*(-P*b - R*T(ii))^3/P^3)^2)/2 - 27*a*b/(2*P) - 9*a*(-P*b - R*T(ii))/(2*P^2) + (-P*b - R*T(ii))^3/P^3)^(1/3)/3 - (-P*b - R*T(ii))/(3*P);
  V3(ii) = -(-3*a/P + (-P*b - R*T(ii))^2/P^2)/(3*(-1/2 + sqrt(3)*I/2)*(sqrt(-4*(-3*a/P + (-P*b - R*T(ii))^2/P^2)^3 + (-27*a*b/P - 9*a*(-P*b - R*T(ii))/P^2 + 2*(-P*b - R*T(ii))^3/P^3)^2)/2 - 27*a*b/(2*P) - 9*a*(-P*b - R*T(ii))/(2*P^2) + (-P*b - R*T(ii))^3/P^3)^(1/3)) - (-1/2 + sqrt(3)*I/2)*(sqrt(-4*(-3*a/P + (-P*b - R*T(ii))^2/P^2)^3 + (-27*a*b/P - 9*a*(-P*b - R*T(ii))/P^2 + 2*(-P*b - R*T(ii))^3/P^3)^2)/2 - 27*a*b/(2*P) - 9*a*(-P*b - R*T(ii))/(2*P^2) + (-P*b - R*T(ii))^3/P^3)^(1/3)/3 - (-P*b - R*T(ii))/(3*P);
end

rho1 = M./abs(V1);
rho2 = M./abs(V3);


plot(T,rho1,T,rho2)
grid on;
grid minor