# lineplot.py
import numpy as np
import pylab as plt

def fermi_dirac(e,mu,beta):
	f=1/(np.exp(beta*(e-mu))+1)
	return f

# Make an x vector array 
x =np.linspace(0,2,100)

# Make an array of y values for each x value
# beta=39.0 implies kT=26 meV  (average KE per particle of ideal gas at 300 K)

y=fermi_dirac(x,1.0,11.6)
w=fermi_dirac(x,1.0,38.7)
z=fermi_dirac(x,1.0,387.0)


#use pylab to plot x and y
plot1 = plt.plot(x, y,'r',label=r'$\beta=11.6$ eV$^{-1}$ (T= 1000K)')
plot2 = plt.plot(x, w,'g',label=r'$\beta=38.7$ eV$^{-1}$ (T=300 K)')
plot3 = plt.plot(x, z,'b',label=r'$\beta= 387$ eV$^{-1}$(T=30K)')
# give plot a title

plt.title('Fermi Dirac Function $\mu$= 1eV')

#make axis labels

plt.xlabel(' Energy (eV)')
plt.ylabel('f(E)')

# set axis limits

#plt.xlim(0.0,7.0)
plt.ylim(0.0,2.0)

plt.legend(loc='upper left')
# show the plot on the screen


plt.show()