import numpy
import pylab

#------------ Parameter definitons ---------------#

# Addition energies, already computed
delta_e = numpy.array([-15.73,-9.15,2.34,8.39])

# Workfunction for the assumed electrodes
w = 5.28

# Gate coupling constant
gate_coupling = 1.0

# Gate bias interval
v_g_interval = numpy.linspace(-15,15,151)

# Source-drain bias interval
v_sd_interval = numpy.linspace(-30,30,301)

#---------- End of parameter definitons ----------#
# test                                           test                         test                test               test                 test
# Calculate the number of charge states in the bias window
def conductionChannels(v_g,v_sd):
    return numpy.sum( abs( delta_e+w+gate_coupling*v_g) <= abs(v_sd/2) )

# Generate the mesh points of the contour plot
X, Y = numpy.meshgrid(v_g_interval,v_sd_interval)

# Evaluate the number of charge states for each mesh point
Z = [ conductionChannels(X[i,j],Y[i,j])
     for i in range(numpy.shape(X)[0])
          for j in range(numpy.shape(X)[1])]

Z = numpy.array(Z).reshape(numpy.shape(X))

# Make the plot
pylab.contour(X,Y,Z)
pylab.contourf(X,Y,Z)
pylab.xlabel("Gate voltage (Volt)")
pylab.ylabel("Source-Drain bias (Volt)")
pylab.show()