%filename Prob_1.m 
clear all
format short e
global press a b R T % make these parameters available to waalsvol.m
%set the constnats
Pcrit=111.3;	% in atm
Tcrit=405.5;	% in Kelvin
R=0.08206;		% in atm.liter/g-mol.K
T=450;		    % K

% the different values of pressure are stored in a single vector
Preduced=[0.503144 1 2 4 10 20];
a=27/64*R^2*Tcrit^2/Pcrit; 
b=R*Tcrit/(8*Pcrit);

% each pass of the loop varies the pressure and the volume is calculated
for j=1:6
  press=Pcrit*Preduced(j);
  volguess=R*T/press;
  % Use fzero ( or fsolve) to calculate volume, the third argument of fzero is t$
  vol= fzero('waalsvol',volguess);
  z=press*vol/(R*T);
  result(j,1)=Preduced(j);
  result(j,2)=vol;
  result(j,3)= press*vol/(R*T);
end
% end of calculation

disp('   Preduced     Molar Vol      Zfactor ')
disp(result)
plot(result(:,1),result(:,3),'r')
title('Compressibility factor vs Reduced pressure')
xlabel('Reduced pressure')
ylabel('Compressibiliity factor')