#!/usr/bin/env python

from pylab import *

T = 300 # K 
kB = 0.02585 / 300 # eV/K
EG = 1.170 - 4.730e-4 * T ** 2 / (T + 636) # pages 69,70
mnr = 1.028 + 6.11e-4 * T - 3.09e-7 * T ** 2 # mnstar / m_e -- page 55
mpr = 0.610 + 7.83e-4 * T - 4.46e-7 * T ** 2 # mpstar / m_e -- page 55

# What do we need to plot?
#   e distribution function, edf = gc(E)*f(E)      -- for E >= Ec
#   h distribution function, hdf = gv(E)*(1-f(E))  -- for E <= Ec
# Define gcp(x) = gc(x+Ec).
# Then,
#   gc(E)*f(E) = gcp(x)*f(x+Ec)
#   where x = E - Ec >= 0 is the electron excitation energy.
# Likewise, define gvp(x) = gv(Ev-x)
#   gv(E)*(1-f(E)) = gvp(x)*(1-f(Ev-x))
#   where x = Ev - E >= 0 is the hole excitation energy.

# x is in unit of eV -- it is positive and used both for e and h.
# We are taking it to go to 7 kB T.
x = linspace (0, 5, 100) * kB * T

# Check that they are reasonable.
print 'mnr =', mnr
print 'mpr =', mpr

# These are from Eqs. T2.6a,b with x = E-Ec or Ev-E.
# If x is in unit of eV, then gcp, gvp are in unit of
#     eV ** 0.5  * eV ** 1.5 / eV ** 3 / Angstrom ** 3
#   = 1 / eV / Angstrom ** 3
gcp = sqrt(2. * x) * (mnr * 0.511e6) ** 1.5 / pi ** 2 / 1973 ** 3
gvp = sqrt(2. * x) * (mpr * 0.511e6) ** 1.5 / pi ** 2 / 1973 ** 3

# As given in the problem.
Ev = 0.
Ec = EG
EF = 0.75 * EG

# Multiply by f(x+Ec) and 1-f(Ev-x) respectively, as shown above.
edf = gcp / (exp ((x + Ec - EF) / kB / T) + 1)
hdf = gvp / (exp ((x - Ev + EF) / kB / T) + 1)

# Finally, convert the unit to 1 / eV / cm ** 3
edf *= 1.e24
hdf *= 1.e24

# Check the orders of magnitude.
print 'max (edf) / max (hdf) = %g' % (max (edf) / max (hdf))

# However, the hole counts are very low -- blow them up.
# Based on what we learnd just now.
hdf *= 5.e9

# All hard compuations now done -- plot them!
plot (x / kB / T, edf, label = 'Electron')
plot (x / kB / T, hdf, label = r'Hole ($\times$  $5\times10^9$)')
grid ()

# matplotlib comes with its own LaTeX parser -- cool!
xlabel (r'Carrier Excitation Energy (in unit of kT)')
ylabel ('Carrier Distribution (1/(eV cm$^3$))')
title ('Carrier Distribution Functions at T = 300 K')
legend ()

# Sort of equivalent to "print" of matlab.  "print" is a python keyword.
# Comment out one of these lines to get a print-out in batch mode.
#savefig ('E2p3.svg')
#savefig ('E2p3.eps')
savefig ('P2p9.png')