#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from scipy.special import gamma
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

k = 1.380649e-23 # J/K
e = 1.602176634e-19 # C
eV = e # J
me = 9.10938356e-31 # kg
h = 6.62607015e-34 # J.s

def calc_A(mass, mobility, a, kT):
    # calculate A based on density-independent parameters.
    # both conductivity and carrier density depend exponentially on
    # chemical potential; their ratio however,
    #   sigma / n = mobility * e
    # has no dependence on chemical potential. So:
    #   A (kT)^a Gamma(a+1) / N = mobility*e
    # where N is concentration coefficient:
    #   N = 2 * (2*pi*mass*kT/h^2)^1.5
    return 2*e*mobility * (2*np.pi*mass)**1.5/(h**3 * gamma(a+1)) * kT**(1.5-a)

T = 300 # K
kT = k*T

# material params for Silicon at 300 K
# band energies
EV = 0*eV
EC = 1.124*eV
EF = np.arange(EV + 1*kT, EC - 1*kT, 0.001*eV)
# effective masses for density of states
m_C = 1.09*me
m_V = 1.15*me
# concentration coeffs (m^-3)
conc_C = 2 * (2*np.pi*m_C*kT/h**2)**1.5
conc_V = 2 * (2*np.pi*m_V*kT/h**2)**1.5
# mobilities
mobility_C = 0.140 # m^2/V/s
mobility_V = 0.045 # m^2/V/s
# calculate conductivity prefactor given known mobility
B_C = mobility_C * conc_C * e
B_V = mobility_V * conc_V * e
# scattering mechanism (acoustic phonon = 1.0)
a_C = 1.
a_V = 1.

# Calculate derived functions
cond_C = B_C * np.exp((EF - EC)/kT)
cond_V = B_V * np.exp((EV - EF)/kT)
n_C = conc_C * np.exp((EF - EC)/kT)
n_V = conc_V * np.exp((EV - EF)/kT)
seeb_C = (-k/e) * (1 + a_C + (EC - EF)/kT)
seeb_V = (+k/e) * (1 + a_V + (EF - EV)/kT)

cond = cond_C + cond_V
seeb = (cond_C * seeb_C + cond_V * seeb_V) / (cond_C + cond_V)


# plotting stuff
# data too near to band edges is bad, make it dashed.
leadin  = (EF < EV + 4*kT)
leadout = (EF > EC - 4*kT)
main = ~(leadin | leadout)

midgap = EF[np.argmax(seeb < 0)] # zero crossing point of seebeck
# midgap = 0.5 * (EV + EC) # halfway

plt.close('all')
fig = plt.figure()
axl = plt.axes()
fig.set_size_inches(4,3)
plt.xlim(EV - 3*kT, EC + 3*kT)
plt.xticks([EV, EC],
           [r"$E_V$", r"$E_C$"])
plt.subplots_adjust(0.15,0.17,0.82,0.98)

axl.plot(EF[main], seeb[main]*1000, color='k', lw=1.5)
axl.plot(EF[leadout], seeb[leadout]*1000, color='k', ls='dashed', lw=1)
axl.plot(EF[leadin], seeb[leadin]*1000, color='k', ls='dashed', lw=1)
axl.set_xlabel(r"Fermi level $\mu$")
axl.set_ylabel(r"Seebeck coefficient $S$ (mV/K)")
#plt.axvline(EV)
#plt.axvline(EC)
#plt.axvline(cross)
anty = -1.99
axl.annotate('',
             xy=(midgap-2*kT, anty), xytext=(midgap+2*kT, anty), xycoords='data', textcoords='data',
             arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0, linewidth=0.5))
axl.annotate(r'$4kT$',
             verticalalignment='center',
             xy=(midgap+2*kT, anty), xytext=(1, 0), xycoords='data', textcoords='offset points')
axl.set_ylim(-2.25, 2.05)
axl.set_yticks([-2,-1,0,1,2])

# draw conductivity content with right axis
axr = axl.twinx()
rcolor='tab:blue'
axr.semilogy(EF[main], cond[main], color=rcolor, lw=1.5, alpha=0.7)
axr.semilogy(EF[leadout], cond[leadout], color=rcolor, ls='dashed', lw=1, alpha=0.7)
axr.semilogy(EF[leadin], cond[leadin], color=rcolor, ls='dashed', lw=1, alpha=0.7)
axr.set_ylabel(r'Conductivity $\sigma$ (S/m)', color=rcolor)  # we already handled the x-label with ax1
axr.tick_params(axis='y', labelcolor=rcolor)
axr.set_ylim(0.1e-4,2e6)
axr.set_yticks([1e-4,1e-2,1,1e2,1e4])

fig.savefig('seeb.svg')