%
%   EQOFST_1.M    Plot equation of state for ammonia over selected range
%                       using van der Waal's equation
%
%   This demo evaluates and plots the equation of state for ammonia over a
%   range of pressures and temperatures.  In this case, the ideal gas law is 
%   compared to van der Waal's equation, which is given by 
%       (P + a/v^2)(v - b) = RT
%   where a and b are constants for a particular gas.  For ammonia, given some
%   value of P and T, the goal is to find the value of molar volume (v in liters/
%   gmole) that satisfies the above equation.  Note that an anonymous function is 
%   used for the implicit form of van der Waal's equation, f(v) = 0
%
%   The compressibility factor, Z = Pv/RT, can be used to compare how close van
%   der Waal's equation is to the ideal gas law.  Note that, for the ideal gas
%   law, Z is always unity, but other (more accurate) forms for the equation of 
%   state will deviate somewhat from unity.
%
%   File prepared by J. R. White, UMass-Lowell (last update:  Nov. 2017)
%

      clear all,   close all,   nfig = 0;
%
%   setup range of variables
      press = [1 3 5];               % pressure (atm)
      temp = linspace(250,400,16)';  % temperature (K)
      R = .08205;                    % universal gas constant (liters-atm/gmole-K)
      Np = length(press);   Nt = length(temp);
%
%   evaluate ideal gas equation (Pv = RT)
      v1 = zeros(Nt,Np);
      for i = 1:Np
        v1(:,i) = R*temp/press(i);
      end
%
%   evaluate van der Waals equation
%   (need a little extra work here because of implicit form)
      v2 = zeros(Nt,Np);   Z = zeros(Nt,Np);
      a = 4.19;  b = .0372;   % constants for ammonia
      for i = 1:Np
        P = press(i);
          for j = 1:Nt
            T = temp(j);
            fv = @(v) (P + a/(v*v))*(v - b) - R*T; % anonymous function used in fzero
            v2(j,i) = fzero(fv,v1(j,i));
            Z(j,i) = P*v2(j,i)/(R*T);
          end
      end
%
%  plot molal volume  -  ideal (dash), van der Waals (solid)
     nfig = nfig+1;   figure(nfig)
     plot(temp,v1,'--',temp,v2,'-','LineWidth',2)
     title('Eqofst\_1:  Equations of State for Ammonia')
     xlabel('Temperature (K)'),ylabel('Molal Volume (liters/gmole)'),grid
     for i = 1:Np,   gtext(['P = ',num2str(press(i)),' atm']),   end
     v = axis;  xt = temp(1)+10;  yt = v(3) + 0.85*(v(4)-v(3));
     txt1 = ['Ideal Gas (dashed line)  '; 'van der Waal (solid line)'];
     text(xt,yt,txt1)
%
%  plot compressibility factor  
     nfig = nfig+1;   figure(nfig)
     plot(temp,Z(:,1),'r-',temp,Z(:,2),'b--',temp,Z(:,3),'g:','LineWidth',2)
     title('Eqofst\_1:  Compressibility Factor for Ammonia (Van der Waals)')
     xlabel('Temperature (K)'),ylabel('Compressibility'),grid
     for i = 1: Np,  txt2(i) = {['P = ',num2str(press(i)),' atm']};  end
     legend(txt2,'Location','SouthEast')
%
%   end of problem
