import numpy as np
import matplotlib.pyplot as plt

plt.rc('text', usetex=True)

def save_twice(name):
    plt.tight_layout()
    plt.savefig('../web/figs/' + name, dpi=600, format='png')
    plt.savefig(name + '.pdf')
    plt.savefig(name + '.svg')
    plt.figure(figsize=(4,3))

b = 1
a = 1
N = 1

n = np.linspace(1e-5, 0.9*N/b, 100000)
V = N/n

fig = plt.figure(figsize=(4,3))

#plt.plot([0,n.max()], [0,0], 'k-')

kT = 0.26*a/b
p = N*kT/(V-N*b) - N**2/V**2*a
plt.plot(n , p, label=r'$kT\frac{b}{a} = %g$' % (kT*b/a))
p0 = 0.0215
okay_points = np.abs(p-p0) < 0.00001
plt.plot([0,n.max()], [p0,p0], 'k:')
plt.plot(n[okay_points], p[okay_points], 'k.')

kT = 0.28*a/b
p = N*kT/(V-N*b) - N**2/V**2*a
plt.plot(n , p, label=r'$kT\frac{b}{a} = %g$' % (kT*b/a))

plt.xlim(0, 0.7)
plt.ylim(min(0, p.min()), 0.07)

kT = 0.3*a/b
p = N*kT/(V-N*b) - N**2/V**2*a
plt.plot(n , p, label=r'$kT\frac{b}{a} = %g$' % (kT*b/a))

plt.ylabel(r"$p\frac{b^2}{a}$")
plt.xlabel(r"$\eta\equiv \frac{N}{V}b$")

plt.legend(loc='best')

save_twice('vdw-pressure')

fig = plt.figure(figsize=(4,3))

kT = 0.26*a/b
p = N*kT/(V-N*b) - N**2/V**2*a
plt.plot(V , p, label=r'$kT\frac{b}{a} = %g$' % (kT*b/a))
okay_points = np.abs(p-p0) < 0.00001
plt.plot([0,V.max()], [p0,p0], 'k:')
plt.plot(V[okay_points], p[okay_points], 'k.')

kT = 0.28*a/b
p = N*kT/(V-N*b) - N**2/V**2*a
plt.plot(V , p, label=r'$kT\frac{b}{a} = %g$' % (kT*b/a))

plt.xlim(0, 20*b)
plt.ylim(min(0, p.min()), .06)

kT = 0.3*a/b
p = N*kT/(V-N*b) - N**2/V**2*a
plt.plot(V , p, label=r'$kT\frac{b}{a} = %g$' % (kT*b/a))

plt.ylabel(r"$p\frac{b^2}{a}$")
plt.xlabel(r"$\frac{V}{Nb}$")

plt.legend(loc='best')

save_twice('vdw-pressure-volume')


fig = plt.figure(figsize=(4,3))

kT = 0.26*a/b
nQ = 5
F = -N*kT*(np.log(nQ*(V-N*b)/N) + 1) - N**2/V*a
plt.plot(V , F, label=r'$kT\frac{b}{a} = %g$' % (kT*b/a))
plt.plot(V[okay_points], F[okay_points], 'k.')

V0 = V[okay_points][0]
V1 = V[okay_points][-1]
F0 = F[okay_points][0]
F1 = F[okay_points][-1]
slope = (F1 - F0)/(V1 - V0)
plt.plot(V, F1 + slope*(V-V1), 'k-')

plt.xlim(0, 13*b)
plt.ylim(-1.4, -1.1)

plt.ylabel(r"$F$")
plt.xlabel(r"$\frac{V}{Nb}$")

plt.legend(loc='best')

save_twice('vdw-free-volume')


fig = plt.figure(figsize=(4,3))

kT = 0.26*a/b
nQ = 5
F = -N*kT*(np.log(nQ*(V-N*b)/N) + 1) - N**2/V*a
p = N*kT/(V-N*b) - N**2/V**2*a
G = F + p*V
plt.plot(p, G, label=r'$kT\frac{b}{a} = %g$' % (kT*b/a))
plt.plot(p[okay_points], G[okay_points], 'k.')

plt.xlim(0.005, .03)
plt.ylim(-1.18, -1.1)

plt.ylabel(r"$G$")
plt.xlabel(r"$p$")

plt.legend(loc='best')

save_twice('vdw-gibbs')