import numpy as np
import matplotlib.pyplot as plt

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

E = np.arange(0, 1.2, 0.0001)
mu = 1
f0 = 1/(np.exp(10000*(E-mu))+1)
f = 1/(np.exp(100*(E-mu))+1)

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


plt.fill_between(E[E<1], f0[E<1], f[E<1], color=[1,.7,.7])
plt.fill_between(E[E>1], f0[E>1], f[E>1], color=[.7,.7,1])

plt.plot(E, f0, label='$T=0$')
plt.plot(E, f, label='$T>0$')
plt.xlabel(r'$\varepsilon$')
plt.ylabel(r'$f(\varepsilon)$')
plt.xlim(0,max(E))
plt.ylim(0, 1.1)

plt.xticks((1,))
fig.axes[0].set_xticklabels((r'$\varepsilon_F$',))
plt.yticks((1,))
fig.axes[0].set_yticklabels(('$1$',))
plt.legend(loc='best')

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))

save_twice('fermi-dirac-distribution')

D = np.sqrt(E)

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

plt.fill_between(E[E<1], D[E<1]*f0[E<1], D[E<1]*f[E<1], color=[1,.7,.7])
plt.fill_between(E[E>1], D[E>1]*f0[E>1], D[E>1]*f[E>1], color=[.7,.7,1])

plt.plot(E, f0*D, label='$T=0$')
plt.plot(E, f*D, label='$T>0$')
plt.xlabel(r'$\varepsilon$')
plt.ylabel(r'$D(\varepsilon)f(\varepsilon)$')
plt.xlim(0,max(E))
plt.ylim(0, 1.1)

plt.xticks((1,))
fig.axes[0].set_xticklabels((r'$\varepsilon_F$',))
plt.yticks([])
fig.axes[0].set_yticklabels([])
plt.legend(loc='best')

save_twice('fermi-gas-product')


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

plt.plot(E, D)
plt.xlabel(r'$\varepsilon$')
plt.ylabel(r'$D(\varepsilon)$')
plt.xlim(0,max(E))
plt.ylim(0, 1.1)

plt.xticks([])
plt.yticks([])
fig.axes[0].set_yticklabels([])
plt.legend(loc='best')

save_twice('gas-density-of-states')