#! /usr/bin/env python

# This will generate the conductance versus surface potential curves
# for undoped Germanium, 15 Ohm-cm n-Germanium, 4 k-cm FZ Si, 
# or 70 Ohm-cm n-Si
# Requires python, numpy, and scipy, which can be obtained from 
# http://www.scipy.org

from math import *
from numpy import *
import csv
import sys

#print sys.argv[1]

# set up the array of sp values
v = linspace( -60.,30.,900 )

while True:
    try:
        print "0 - Conductance Curve, 1 - Conversion to Potential"
        purp = int(raw_input('-->'))
        break
    except ValueError:
        print "Try again..."

while True:
    try:
        print "0 - Undoped Ge, 1 - 15 Ohm-cm n-Ge, 2 - 4k-cm FZ Si, 3 - 70 Ohm-cm Si" 
        dope = int(raw_input('-->'))
        break
    except ValueError:
        print "Try again..."

if dope <= 1:
# parameters for Ge
    kappa = 16
    e0 = 8.85e-014 # F/cm
    q = 1.60e-19 #C
    ni = 2.40e+13 #cm-3
    k = 8.61e-05 #eV/K
    L = 6.8e-05 #cm debye, intrinsic
    T = 300 # K
    mun = 3800 #
    mup = 1820 #
else:
# parameters for Si

    kappa = 11.9
    e0 = 8.85e-014 # F/cm
    q = 1.60e-19 #C
    ni = 1.45e+10 #cm-3
    k = 8.61e-05 #eV/K
    L = 2.4e-03 #cm debye, intrinsic
    T = 300 # K
    mun = 1500 # cm2/V-s
    mup = 450 #


if dope == 0:
# Ge-Intrinsic
# gmind is for minority carriers under depletion/inversio
# gmaja is for majority carriers under accumulation, etc.


    def gp(x):
        return exp(0.5*abs(x))-1

    def gn(x):
        return exp(-0.5*abs(x))-1

    gmind = 2*L*ni*gp(v)
    gmina = 2*L*ni*gn(v)

# calling electrons majority, even for intrinsic

    gmaja = gmind
    gmajd = gmina


# These values use:
# nb = impurity concentration (Hall or fig.22 Ch1 Sze)

elif dope == 1:
# Ge-15 Ohm-cm
    nb = 1.0e+14
    pb = (ni**2)/nb
    ub = log(nb/ni)
    Le = L*(((2*ni)/(nb+pb))**0.5) #effective Debye length 

elif dope == 2:
# Si-4 k-cm
    nb = 4.5e+11
    pb = (ni**2)/nb
    ub = log(nb/ni)
    Le = L*(((2*ni)/(nb+pb))**0.5)

elif dope == 3:
# Si-70 Ohm-cm
    nb = 6.6e+13
    pb = (ni**2)/nb
    ub = log(nb/ni)
    Le = L*(((2*ni)/(nb+pb))**0.5)

# now calculate for extrinsic
if dope == 1 or dope == 2 or dope == 3:
    def fuv1(ub,v1):
        return cosh(ub+v1)/cosh(ub)
    def fuv2(ub,v1):
        return 2.*fuv1(ub,v1)-2.*v1*tanh(ub)-2
    def fuv3(ub,v1):
        return fuv2(ub,v1)**0.5
# majority carriers, accumulation
    def integrand_maja(v1):
        return (exp(v1)-1.)/fuv3(abs(ub),v1)
# majority carriers, depletion
    def integrand_majd(v1):
        return (exp(-v1)-1.)/fuv3(-abs(ub),v1)
# minority carriers, accumulation
    def integrand_mina(v1):
        return exp(-2.*abs(ub))*(exp(-v1)-1.)/fuv3(abs(ub),v1)
# minority carriers, depletion
    def integrand_mind(v1):
        return exp(-2.*abs(ub))*(exp(v1)-1.)/fuv3(-abs(ub),v1)
# integration
    from scipy.integrate import quad

    def gna(ub,vs):
        return quad(integrand_maja,0.,abs(vs),args=())[0]
    def gnd(ub,vs):
        return quad(integrand_majd,0.,abs(vs),args=())[0]
    def gpa(ub,vs):
        return quad(integrand_mina,0.,abs(vs),args=())[0]
    def gpd(ub,vs):
        return quad(integrand_mind,0.,abs(vs),args=())[0]

    gmaja = zeros((900),dtype=float)
    gmina = zeros((900),dtype=float)
    gmind = zeros((900),dtype=float)
    gmajd = zeros((900),dtype=float)
    it=0
    while it < 900:
        gmaja[it] += Le*nb*gna(ub,abs(v[it]))
        gmina[it] += Le*nb*gpa(ub,abs(v[it]))
        gmajd[it] += Le*nb*gnd(ub,abs(v[it]))
        gmind[it] += Le*nb*gpd(ub,abs(v[it]))
        it = it + 1



# this is more generalized
# array 'sig' has two axes, 600 is for each surface potential
# 2 is for (x-sp,y-sc)  
sig = zeros((900,2), dtype=float)
# only add to depleted/inverted
sig[0:599,1] += mun*gmajd[0:599]
sig[0:599,1] += mup*gmind[0:599]
# only add to accumulated
sig[600:899,1] += mun*gmaja[600:899]
sig[600:899,1] += mup*gmina[600:899]
# multiply by q and we are done
sig[:,1] *= q
sig[:,0] += v

if purp == 0:
    fest = csv.writer(open('/home/dk246/Thesis/SPmap_pub/curveout-dummy.csv','w'),delimiter=' ')
    print fest
    it=0
    while it < 899:
        fest.writerow([sig[it,0],sig[it,1]])
        it = it + 1
####################################################################

if purp == 1:

# find the minimum
    it = 0
    cmin = 1.0e-03 # an unlikely value
    while it < 899:
        if sig[it,1] < cmin:
            cmin = sig[it,1]
            pmin = sig[it,0]
        it = it + 1  

# now the data actually gets read
# The file to read

    filein = '/home/dk246/Thesis/SPmap/SPmap_in/%s.csv' % (sys.argv[1])
    fileout = '/home/dk246/Thesis/SPmap/SPmap_out/%s.csv' % (sys.argv[1])
    f = open(filein,'r')
    af=f.readlines()
    l=0 # line number
    data_out = zeros((len(af),2), dtype=float)
    apot_old = -50. # so that we choose the correct arm of curve
#                     it needs to be changed for the oxide
    a_old = 1.0e-3
    while l < len(af): # len(af) is from readlines
        ax = af[l]
        it=0
        while it < 20:  # 20 is just a nice size
            it += 1
            if ax[it] == ',':
                break
        a = float(ax[0:it]) # a is the conductacnce
# the data is conductance change relative to minimum, but calculations are relaive to flat-band
        a = a + cmin
        b = float(ax[it+1:-2]) # b is the SRV
        it = 0
        ahold1 = 1.0e-03 # some unlikely residual
        ahold2 = 1.0e-03 # some unlikely residual
        apot_out = 0.
        while it < 899:
            adiff = a - sig[it,1]
            if abs(adiff) < abs(ahold1) and sig[it,0] < pmin:
                apot1 = sig[it,0]
                ahold1 = adiff
            elif abs(adiff) < abs(ahold2) and sig[it,0] > pmin:
                apot2 = sig[it,0]
                ahold2 = adiff
            it = it + 1

# two unambiguous cases
# conductance is decreasing, and old potential is closer to negative choice
        if a < a_old and abs(apot_old - apot1) < abs(apot_old - apot2): 
            apot_out = apot1
            apot_old = apot1
            a_old = a
# conductance has increased, and old potential is closer to positive choice
        elif a >= a_old and abs(apot_old - apot1) >= abs(apot_old - apot2):    
            apot_out = apot2
            apot_old = apot2
            a_old = a
# ambiguous cases
# the coductance has increased, but last potential is closer to negative choice
# choose negative potential, but store positive potential and don't overwrite the old conductance.  
# If it is not a jitter, the next one will be increased too 
# and stored potential will be closer to positive choice next time
# If it is jitter, the next conductance will not be increased and will be dealt with below
        elif a >= a_old and abs(apot_old - apot1) < abs(apot_old - apot2):    
            if abs(abs(apot_old - apot1) - abs(apot_old - apot2)) < 3.0: # about 75 mV
                apot_out = apot2 # it is a toss-up, just proceed
                apot_old = apot2
            else:
                apot_out = apot1
                apot_old = apot2
# The conductance has decreased, but last potential is closer to poitive choice
# choose positive potential but store negative and don't overwrite conductance
# If it is not jitter, next conductance will decrease too, and stored potential 
# will be closer to negative next time
# If it is jitter, next conductance will increase but stored potential will be 
# closer to negative, so it will go to case above for one value, then get back on track 
        else:
            if abs(abs(apot_old - apot1) - abs(apot_old - apot2)) < 3.0:
                apot_out = apot1
                apot_old = apot1
            else:
                apot_out = apot2
                apot_old = apot1


# so that output is in volts
        apot_out *= k*T
        data_out[l,0] = apot_out
        data_out[l,1] = b
        l = l + 1
    it = 0
    f2 = csv.writer(open(fileout,'w'),delimiter=' ')
    while it < len(af):
        f2.writerow([data_out[it,0],data_out[it,1]])
        it = it + 1