%     van der Waals Adsorption Isotherm
%     a2 cohesive parameter, J m^2
%     b2 coarea, m^2
%     A  Hamaker constant, J
%     hi interaction distance, m
%     k  Boltzmann's constant, J/K
%     pr p/P(T)
%     p  partial pressure, mm Hg
%     PT vapor pressure at T, mm Hg
%     theta fractional coverage, b2/a
%     T  temperature, K
%     Tc2 2D critical temperature, K
%     v  molecular volume m^-3
%     H  d(theta)/d(pr), dimensionless Henry's law constant
      clear all
      figure(1)
      clf reset
      figure(2)
      clf reset
      figure(3)
      clf reset
      load vdwdata.dat
      A=vdwdata(1);
      k=vdwdata(2);
      v=vdwdata(3);
      hi=vdwdata(4);
      a2=vdwdata(5);
      b2=vdwdata(6);
      Tc2=(8/27)*a2/(k*b2)
      niso=input('Number of Isotherms = ')
      for n=1:niso
      Tcel=input('T, deg C = ');
      T=Tcel+273.16;
      PT=input('Vapor Pressure @ T, mm Hg =');
      C1=(1/(6*pi))*(A/(k*T))*v*hi^(-3);
      H=exp(-C1);
      ['Dimensionless Henrys law constant = ', num2str(H)]
      theta=(0.01:0.01:0.98)';
      pr=zeros(size(theta));
      pr(:)=(theta(:)./(1-theta(:))).*exp(C1-(27/4)*(Tc2/T)*theta(:) ...
            +(theta(:)./(1-theta(:))));
      figure(1)
      semilogx(pr,theta)
      axis([1.e-6 1 0 1])
      xlabel('p/P(T)')
      ylabel('Theta')
      if niso==1 title(['T = ',num2str(T),' deg. K,  Vap. Press. = ',num2str(PT),...
            ' mm Hg'])
      end
      hold on
      figure(2)
      p=PT.*pr;
      semilogx(p,theta)
      axis([1.e-6 1.e3 0 1])
      xlabel('p, mm Hg')
      ylabel('Theta')
      if niso==1 title(['T = ',num2str(T),' deg. K,  Vap. Press. = ',num2str(PT),...
            ' mm Hg'])
      end
      hold on
      figure(3)
      plot(p,theta)
      axis([0 40 0 1])
      xlabel('p, mm Hg')
      ylabel('Theta')
      if niso==1 title(['T = ',num2str(T),' deg. K,  Vap. Press. = ',num2str(PT),...
            ' mm Hg'])
      end
      hold on
      end  % End of do loop
