import numpy as np
import matplotlib.pyplot as plt


# The density of states (dos) function to plot
def dos(E,D):
    if np.abs(E)>D:
        return np.abs(E)/np.sqrt(E**2-D**2)
    else:
        return 0


# Create lists of 'x' and 'y' values to plot
Elist = np.arange(-4,4,0.01)
doslist = []

for Eval in Elist:
    doslist.append(dos(Eval,1))



# Some plot variables
Elabel = 'E/$\Delta$'
Doslabel = 'N'

E1 = -3
E2 = 3

Dos1 = 0
Dos2 = 3

outfile = 'basic_plot1.png'



# Create the plot
fig = plt.figure(figsize=(8,5))
ax1 = fig.add_subplot(1,1,1)
ax1.plot(Elist,doslist, label=1, linewidth=1, linestyle="-", color='blue')


# Set axis parameters
ax1.legend(frameon=False,title='$\Delta$',fontsize='large',loc=1)
ax1.set_xlabel(Elabel)
ax1.set_ylabel(Doslabel)
ax1.set_xlim(E1,E2);
ax1.set_ylim(Dos1,Dos2);
    

# Save the plot
plt.savefig(outfile, format="png", bbox_inches='tight')        
plt.show()