import csv
import numpy as np
from numpy import sin, cos, tan, pi, exp
import sys

# Brent's method
# Algorism and virtual algorism are referred from
#    https://ja.wikipedia.org/wiki/%E3%83%96%E3%83%AC%E3%83%B3%E3%83%88%E6%B3%95

# physical parameeters
e  = 1.60218e-19      # C";
kB = 1.380658e-23     # JK<sup>-1</sup>";

T = 300.0
ekBT = e / kB / T

# semiconductor parameters
Ev = 0.0
Ec = 1.0
Nv = 1.2e19
Nc = 2.1e18
EA = 0.05
NA = 1.0e15
ED = Ec - 0.05
ND = 1.0e16

# parameters
EFmin = Ev
EFmax = Ec
eps   = 1.0e-5
#delta = (EFmax - EFmin) * 1.0e-1
delta = eps
nmaxiter = 100
iprintinterval = 1

# Fermi-Dirac function
def fe(E, EF, T):
    global Nc, Ec, kB
    return 1.0 / (exp((E - EF) * ekBT) + 1.0)

# electron density
def Ne(EF, T):
    global Nc, Ec, kB
    return Nc * exp(-(Ec - EF) * ekBT)

# hole density
def Nh(EF, T):
    global Nv, Ev, kB
    return Nv * exp(-(EF - Ev) * ekBT)

# ionized donor density
def NDp(EF, T):
    global ND, ED, kB
    return ND * (1.0 - fe(ED, EF, T))

# ionized acceptor density
def NAm(EF, T):
    global NA, EA, kB
    return NA * fe(EA, EF, T)


def main():
    global EFmin, EFmax, delta, eps, nmaxiter, iprintinterval

    print("Solution of EF by Brent's method")

    print("")
    [EFa, EFb] = [EFmin, EFmax]
    dQa = Ne(EFa, T) + NAm(EFa, T) - Nh(EFa, T) - NDp(EFa, T)
    dQb = Ne(EFb, T) + NAm(EFb, T) - Nh(EFb, T) - NDp(EFb, T)
    if dQa * dQb > 0.0:
        print("Error: Initial Emin and Emax should be chosen as dQmin * dQmax < 0")
        return 0
    if abs(dQa) < abs(dQb):
        [EFa, EFb] = [EFb, EFa]
        [dQa, dQb] = [dQb, dQa]
    print("  EFmin = {:12.8f}  dQmin = {:12.4g}".format(EFa, dQa))
    print("  EFmax = {:12.8f}  dQmax = {:12.4g}".format(EFb, dQb))

    EFc = EFa
    mflag = 1

    for i in range(nmaxiter):
        if abs(EFa - EFb) < eps:
            EFh = (EFa + EFb) / 2.0
            print("  Success: Convergence reached at EF = {}".format(EFh))
            return 1

        dQc = Ne(EFc, T) + NAm(EFc, T) - Nh(EFc, T) - NDp(EFc, T)
        if abs(dQa - dQc) > 1.0e-10 and abs(dQb - dQc) > 1.0e-10:
#inverse quadratic interpolation
            EFs = EFa * dQb * dQc / (dQa - dQb) / (dQa - dQc) \
                + EFb * dQa * dQc / (dQb - dQa) / (dQb - dQc) \
                + EFc * dQa * dQb / (dQc - dQa) / (dQc - dQb)   
        else:
#secant method
            EFs = EFb - dQb * (EFb - EFa) / (dQb - dQa)

        EF4 = (3.0*EFa + EFb) / 4.0
        if not (EF4 < EFs < EFb or EFb < EFs < EF4):
            mflag = 1
#            print("condition 1")
        elif mflag == 1 and abs(EFs - EFb) >= abs(EFb - EFc) / 2.0:
            mflag = 1
#            print("condition 2")
        elif mflag == 0 and abs(EFs - EFb) >= abs(EFc - EFd) / 2.0:
            mflag = 1
#            print("condition 3")
        elif mflag == 1 and abs(EFb - EFc) < delta:
            mflag = 1
#            print("condition 4")
        elif mflag == 0 and abs(EFc - EFd) < delta:
            mflag = 1
#            print("condition 5")
        else:
            mflag = 0
#            print("condition 6")
        if mflag == 1:
            EFs = (EFa + EFb) / 2.0  # bisection

        dQs = Ne(EFs, T) + NAm(EFs, T) - Nh(EFs, T) - NDp(EFs, T)
        EFd = EFc
        dQd = dQc
        EFc = EFb
#        dQc = dQb

        if dQa * dQs < 0.0:
#            print("condition a: EFa,b,c,s=", EFa, EFb, EFc, EFs)
            EFb = EFs
            dQb = dQs
        else:
#            print("condition b: EFa,b,c,s=", EFa, EFb, EFc, EFs)
            EFa = EFs
            dQa = dQs
        if abs(dQa) < abs(dQb):
            [EFa, EFb] = [EFb, EFa]
            [dQa, dQb] = [dQb, dQa]

        Nes  = Ne(EFs, T)
        NAms = NAm(EFs, T)
        Nhs  = Nh(EFs, T)
        NDps = NDp(EFs, T)
        dQs  = Nes + NAms - Nhs - NDps
        if i % iprintinterval == 0:
            print("  Iter {}: EFa,b = {:12.8f} - {:12.8f}  dQa,b = {:12.4g} - {:12.8g}".format(
                        i, EFa, EFb, dQa, dQb))
            print("     EFs = {:12.8f}  dQs = {:12.4g}".format(i, EFs, dQs))
            print("     Ne={:10.4e}  Nh={:10.4e}  NA-={:10.4e}  ND+={:10.4e}  dQ={:10.4e}".format(
                        Nes, Nhs, NAms, NDps, dQs))

    else:
        print("  Failed: Convergence did not reach")
        return 0


if __name__ == "__main__":
    main()
