Bragg_Williams.m

% Author: Housam Binous

% Pressure-volume isotherms for a Bragg-Williams lattice gas

% National Institute of Applied Sciences and Technology, Tunis, TUNISIA

% Email: [e-mail]

% Reference: T. L. Hill, An Introduction to Statistical Thermodynamics,
% Dover, New York, 1986.

clc

global P1

% Plotting Isotherms. We get result similar to the van der Waals equation.
% The Bragg-Williams approximation pedicts a first-order phase transition.

figure(1)
hold on

for kTovercw=-0.5:0.1:-0.1
    u=1.0001:0.01:10;
    P=kTovercw*log(1 - 1./u) - 1./(2*u.^2);
    h=plot(u,P,'r');
    set(h,'color',rand(3,1))
end
axis([0 10 -0.3 0.3])

% We use the symmetry condition: taul+taug=1 to draw the horizontal stable
% equilibrium path. We can also use the thermodynamic equal-area 
% theorem shown below.

X=fsolve(@zeroBragg,[1.1 5.5])

kTovercw=-0.2;
P1=kTovercw*log(1 - 1./X(1)) - 1./(2*X(1).^2)
P2=kTovercw*log(1 - 1./X(2)) - 1./(2*X(2).^2)

Y=fsolve(@(u) kTovercw*log(1 - 1./u) - 1./(2*u.^2)-P1, 2)
P3=kTovercw*log(1 - 1./Y) - 1./(2*Y.^2)

figure(2)
hold on
u=1.0001:0.01:10;
P=kTovercw*log(1 - 1./u) - 1./(2*u.^2);
h=plot(u,P,'r');
set(h,'color',rand(3,1))
h=plot(X(1),P1,'b.');
set(h,'markersize',20)
h=plot(X(2),P2,'b.');
set(h,'markersize',20)
h=plot(Y,P3,'b.');
set(h,'markersize',20)
axis([0 10 -0.02 0.1])

x=X(1):0.01:Y;
b=kTovercw*log(1 - 1./x) - 1./(2*x.^2);
a=P1*ones(1,122);
[ph,msg]=jbfill(x,a,b,rand(1,3),rand(1,3),0,rand(1,1));

x=Y:0.01:X(2);
b=kTovercw*log(1 - 1./x) - 1./(2*x.^2);
a=P1*ones(1,453);
[ph,msg]=jbfill(x,a,b,rand(1,3),rand(1,3),0,rand(1,1));

% As expected, both shaded areas are equal (A1=A2)

A1=-quadl(@pBragg,X(1),Y)

A2=quadl(@pBragg,Y,X(2))