%Define known constants
R=0.08206;
TC=647.4;
PC=218.3; 
a0=7.700;
a1=0.04594e-2;
a2=0.2521e-5;
a3=-0.8587e-9;
hl0=0;
T0=273.15;
hvap0=(2501.3*18)/1000;
a=0.42747*((R^2)*(TC^(5/2))/PC);
b=0.08664*(((R)*(TC))/PC);
 
%Define inputs
T=373;
P = 1;
if T<0
    disp('T cannot be less than 0 K')
end
if P<=0
    disp('P cannot be less than 0 atm')
end
 
%Use Reidlich Kwong equation to solve for V by using Newtown Raphson Root
%Finding Method
V=10; %Initial guess (place to start while loop)
error=10; %arbitrary initial error
IT=0; %Counts iterations
 
while error>=10^-5 %while loop continues until the difference between V and Vguess is greater than 10^-5
IT=IT+1;
S=((R*T)/(V-b))-(a/(V*(V+b)*sqrt(T)))-P;
Sprime=((-R*T)/(V-b)^2)+(((2*V+b)*a)/((V^2+V*b)^2*sqrt(T)));
Vnew=V-(S/Sprime);
error=abs(Vnew-V);
V=Vnew;
end
 
V;
disp('V (m^3/kg-mol) =')
disp(V)
 
Vs1=V/18;
disp('Vs1 (m^3/kg) =')
disp(Vs1)
 
n=1;
%VI=(n*R*T)/P;
%z=(VI-V)/V
z1=(P*V)/(R*T);
disp('z1 = ')
disp(z1)
 
syms t
q=int(a0+(a1*t)+(a2*t^2)+(a3*t^3),T0,T);
Q=q/(0.23901*1000);
 
Hv1=hvap0+Q+(R*T)*(z1-1-((1.5*a)/(b*R*(T^(1.5))))*log(1+(b/V)));
disp('Hv1 (kJ/mol) =')
disp(Hv1)


