import matplotlib.pyplot as plt
import numpy as np
import scipy.constants as c
e = c.physical_constants["electron volt"][0]
k = c.k # boltzman constant
u = 0 # chemical potential (eV)
E = np.linspace(-.1, .1, 1000) # Energy in eV
T = 10 # temp
a = lambda energy: np.exp(((energy - u) * e) / (c.k * T))
b = lambda energy: np.sqrt(energy) * np.exp(((energy - u) * e) / (c.k * T))
n_b = lambda e: 1 / a(e)
n_fd = lambda e: 1 / (a(e) + 1)
n_be = lambda e: 1 / (a(e) - 1)
n_b_s = lambda e: 1 / b(e)
n_fd_s = lambda e: 1 / (b(e) + 1)
n_be_s = lambda e: 1 / (b(e) - 1)
fig, ax = plt.subplots(1, 1)
ax = plt.gca()
ax.set_xlim([-0.01, 0.01])
ax.set_ylim([0, 1])
ax.plot(E, n_b(E), 'g-', lw=5, alpha=0.6, label='maxwell-boltzmann')
ax.plot(E, n_fd(E), 'r-', lw=5, alpha=0.6, label='fermi-dirac')
ax.plot(E, n_be(E), 'b-', lw=5, alpha=0.6, label='bose-einstein')
plt.legend(loc="lower left")
fig, ax = plt.subplots(1, 1)
ax = plt.gca()
ax.set_xlim([0, 0.01])
ax.set_ylim([0, 1.2])
ax.plot(E, n_b_s(E), 'g-', lw=5, alpha=0.6, label='maxwell-boltzmann')
ax.plot(E, n_fd_s(E), 'r-', lw=5, alpha=0.6, label='fermi-dirac')
ax.plot(E, n_be_s(E), 'b-', lw=5, alpha=0.6, label='bose-einstein')
plt.legend(loc="lower left")