import os
import sys
from math import sqrt, log, exp
import numpy as np
from matplotlib import pyplot as plt
import pandas as pd


from tkpnjunction import tkPNJunction


kB = 1.380649e-23     # JK<sup>-1</sup>";
e0 = 8.854418782e-12; # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
e  = 1.602176634e-19  # C";
pi = 3.14159265358979323846


class tkParams():
    def __init__(self):
        pass


def initialize():
    cparams = tkParams()

    cparams.mode = 'NW'
    cparams.plot = 1

#    cparams.infile  = 'C-V評価 生データ.xlsx'
    cparams.infile  = 'CV_meas.xlsx'
    cparams.outfile = 'CV.xlsx'

    cparams.T   = 300.0 # K
#    r_electrode = 1050e-6  # 850e-6 m
    r_electrode = 1095e-6  # 850e-6 m
    cparams.S   = r_electrode * r_electrode# m2
    cparams.Vbi = 1.0  # V
    cparams.V   = 0.0  # V

# N+ layer
    cparams.dnn   = 7.0e-6    # m
    cparams.NCnn  = 5.0e18    # cm-3
    cparams.NVnn  = 1.0e20    # cm-3
    cparams.NDnn  = 1.0e19    # cm-3
    cparams.epsnn = 10        # e0

# N layer
    cparams.dn   = 3.1e-6    # 3.7e-6 # m
    cparams.NCn  = 5.0e18    # cm-3
    cparams.NVn  = 1.0e20    # cm-3
    cparams.ND   = 1.2e16 #5.0e16    # 2.7e16 # cm-3
    cparams.epsn = 10        # e0

# P layer
    cparams.dp   = 375e-9 # m
    cparams.NCp  = 5.0e18    # cm-3
    cparams.NVp  = 1.0e20    # cm-3
    cparams.NA   = 1.8e18  # 5.0e18 # cm-3
    cparams.epsp = 10     # e0

    cparams.Vmin = -800 # V
    cparams.Vmax = 0    # V
    cparams.nV   = 1001

    cparams.figsize             = [8, 6]
    cparams.fontsize            = 14
    cparams.legend_fontsize     = 12

    return cparams

def update_variables(cparams):
    argv = sys.argv
    nargs = len(argv)
    
    if nargs >= 2:
        cparams.mode = argv[1]
    if cparams.mode == 'NW':
        if nargs >= 3:
            cparams.infile = argv[2]
        if nargs >= 4:
            cparams.outfile = argv[3]
        if nargs >= 5:
            cparams.S    = float(argv[4])
        if nargs >= 6:
            cparams.ND   = float(argv[5])
        if nargs >= 7:
            cparams.NA   = float(argv[6])
        if nargs >= 8:
            cparams.dn   = float(argv[7])
        if nargs >= 9:
            cparams.dp   = float(argv[8])
        if nargs >= 10:
            cparams.plot = int(argv[9])
    else:
        print()
        print(f"Invalid mode [{cparams.mode}]")
        input("Press ENTER to terminate>>")
        exit()

def read_file(infile, cparams):
    if not os.path.isfile(infile):
        return None, None, None

    df = pd.read_excel(infile, header = None)
    v00 = df.iloc[0,0]
    if type(v00) != str or v00 == '':
        labels = df.iloc[1, 1:].tolist()
        data = df.iloc[2:, 1:].values
    else:
        labels = df.iloc[0, 0:].tolist()
        data = df.iloc[1:, 0:].values

    data_df = pd.DataFrame(data, columns = labels)
#    pd.options.display.float_format = '{:.5g}'.format
    print()
    print(f"Read data from [{infile}]")
    labels = data_df.columns.tolist()
    data_list = data_df.values.tolist()
    ndata = data_df.shape[0]
    nskip = max(1, int(ndata / 25))
    print("ndata=", ndata, nskip)
    print(f"{labels[0]:>10} {labels[1]:>10} {labels[1] + '(F/m2)':>10}")
    V_list   = [data_list[i][0] for i in range(ndata)]
    C_list   = [data_list[i][1] for i in range(ndata)]
    Cm2_list = [data_list[i][1] / cparams.S for i in range(ndata)]
#    print(data_list)
    for i in range(0, ndata, nskip):
        print(f"{V_list[i]:10.4g} {C_list[i]:10.4g} {Cm2_list[i]:10.4g}")
        
    return V_list, C_list, Cm2_list


def plot_NW(cparams):
    print()
    print("Calculate N - W relation for n+/n/p junction")
    T = cparams.T
    print(f"T={cparams.T} K")
    print(f"S={cparams.S} m^2")
    print(f"V={cparams.V} V")
    print(f"N+ layer:")
    EFnn = kB * T / e * log(cparams.NDnn / cparams.NCnn)
    print(f"  d  ={cparams.dnn * 1.0e6} um")
    print(f"  NC ={cparams.NCnn} cm^-3")
    print(f"  ND ={cparams.NDnn} cm^-3")
    print(f"  eps={cparams.epsnn} e0")
    print(f"    EF - EC ={-EFnn} eV")
    print(f"N layer:")
    EFn = kB * T / e * log(cparams.ND / cparams.NCn)
    print(f"  d  ={cparams.dn * 1.0e6} um")
    print(f"  NC ={cparams.NCn} cm^-3")
    print(f"  ND ={cparams.ND} cm^-3")
    print(f"  eps={cparams.epsn} e0")
    print(f"    EF - EC ={-EFn} eV")
    print(f"P layer:")
    EFp = kB * T / e * log(cparams.NA / cparams.NVp)
    print(f"  d  ={cparams.dp * 1.0e6} um")
    print(f"  NV ={cparams.NVp} cm^-3")
    print(f"  NA ={cparams.NA} cm^-3")
    print(f"  eps={cparams.epsp} e0")
    print(f"    EF - EV ={-EFp} eV")

    Vbi_nn_n    = EFnn - EFn
    Vbi_n_p     = cparams.Vbi - Vbi_nn_n
    print("Junctions:")
    print(f" Vbi(total)={cparams.Vbi} eV")
    print(f"  Vbi(N+/N)={Vbi_nn_n:.4f} eV")
    print(f"  Vbi(N/P)(to be used) ={Vbi_n_p:.4f} eV")
    W_nn_dep = sqrt(2.0 * cparams.epsnn * e0 * Vbi_nn_n / e / (cparams.NDnn * 1.0e6) )  # m
    Cnn = cparams.epsnn * e0 / W_nn_dep     # F/m2
    W_n_dep  = sqrt(2.0 * cparams.epsn  * e0 * Vbi_n_p  / e / (cparams.ND   * 1.0e6) )  # m
    W_p_dep  = sqrt(2.0 * cparams.epsn  * e0 * Vbi_n_p  / e / (cparams.NA   * 1.0e6) )  # m
    print(f"  W(N+, depletion) at Vbi={Vbi_nn_n:.4g} eV={W_nn_dep * 1.0e9:10.4g} nm")
    print(f"  C(N+, depletion) at Vbi={Vbi_nn_n:.4g} eV={Cnn:10.4g} F/m^2")
    print(f"  W(N, depletion)  at Vbi={Vbi_n_p:.4g} eV={W_n_dep * 1.0e9:10.4g} nm")
    print(f"  W(P, depletion)  at Vbi={Vbi_n_p:.4g} eV={W_p_dep * 1.0e9:10.4g} nm")
    
    print()
    print( "#===============================================================================")
    print(f"# Calculation V range: {cparams.Vmin} - {cparams.Vmax} V, {cparams.nV} points")
    print( "#===============================================================================")

    print()
    print(f"Read data from [{cparams.infile}]")
    V_list, C_meas, Cm2_meas = read_file(cparams.infile, cparams)

    if V_list is not None:
        nV = len(V_list)
        Cm2inv_meas = [1.0 / Cm2_meas[i] / Cm2_meas[i] for i in range(nV)]
        V = V_list
    else:
        nV = cparams.nV
        C_meas = None
        Cm2_meas = None
        Cm2inv_meas = None
        Vstep = (cparams.Vmax - cparams.Vmin) / (nV - 1)
        V = [cparams.Vmin + i * Vstep for i in range(nV)]

    nskip = int(nV / 50)

    Cnmin  = cparams.epsn * e0 / cparams.dn     # F/m2
    V_Cp = []
    Cp   = []
    Cp2inv = []
    print()
    print(f"Extract p-layer Cp using full-depleted n Cp={Cp} F/m2 & Cnn={Cnn} F/m2")
    print(f"{'V (V)':8}\t{'Cmeas (F/m2)':10}\t{'Cp (F/m2)':10}")
    for iV, _V in enumerate(V_list):
        _Cp_inv = 1.0 / Cm2_meas[iV] - 1.0 / Cnmin - 1.0 / Cnn
#        print(iV, _V, Cm2_meas[iV], Cnmin, Cnn, 1.0 / _Cp_inv)
        if _Cp_inv <= 0.0:
            continue
        _Cp = 1.0 / _Cp_inv

        V_Cp.append(_V)
        Cp.append(_Cp)
        Cp2inv.append(1.0 / _Cp / _Cp)
        if iV % nskip == 0:
            print(f"{_V:8.4g}\t{Cm2_meas[iV]:10.4g}\t{_Cp:10.4g}")

    N_Cp = []
    W_Cp = []
    n_Vp = len(V_Cp)
    print()
    print(f"{'V (V)':8}\t{'Cp (F/m2)':10}\t{'Wp (nm)':10}\t{'N (cm-3)':10}")
    for iV, _V in enumerate(V_Cp):
        if iV == 0:
            ip = 1
            im = 0
        elif iV == n_Vp - 1:
            ip = n_Vp - 1
            im = n_Vp - 2
        else:
            ip = iV + 1
            im = iV - 1
        diff = (Cp2inv[ip] - Cp2inv[im]) / (V_Cp[ip] - V_Cp[im])
        _N = -2.0 / e / cparams.epsp / e0 / diff   # m^-3
        _W = cparams.epsp * e0 / Cp[iV]

        N_Cp.append(_N)
        W_Cp.append(_W)
        print(f"{_V:8.4g}\t{Cp[iV]:10.4g}\t{_W:10.4g}\t{_N*1.0e-6:10.4g}")
    
    pnj = tkPNJunction(cparams.ND, cparams.epsn, cparams.dn, cparams.NA, cparams.epsp, cparams.dp)

    C = []
    C2inv = []
    Wn = []
    Wp = []
    W  = []
    Vn = []
    Vp = []
    Vstep = (cparams.Vmax - cparams.Vmin) / (nV - 1)
    nskip = int(cparams.nV / 50)
    Vn_fdep = None 
    for i in range(nV):
        _V = V[i]
#        _Wn, _Wp, _W, _Vn, _Vp, _Vtot, _Cn, _Cp, _C = pnj.cal_depletion_layers(_V, cparams.Vbi)
        _Wn, _Wp, _W, _Vn, _Vp, _Vtot, _Cn, _Cp, _C = pnj.cal_depletion_layers(_V, Vbi_n_p)
        Ctot = 1.0 / (1.0 / Cnn + 1.0 / _Cn + 1.0 / _Cp)
        _C = Ctot

        if Vn_fdep is None and _Wn >= cparams.dn:
            Vn_fdep = _V   

        C.append(_C)
        _C2inv = 1.0 / _C / _C
        C2inv.append(_C2inv)

        Wn.append(_Wn)
        Wp.append(_Wp)
        W.append(_Wn + _Wp)
        Vn.append(_Vn)
        Vp.append(_Vp)
#        if i % nskip == 0:
#            print(f"{_V:10.4g} {_Wn*1e9:10.4g} {_Wp*1e9:10.4g} {_C:10.4g} {_C2inv:10.4g}")

    N = []
    diff_meas = []
    N_meas = []
    W_meas = []
    print()
    print("nV=", nV, nskip)
    if C_meas:
        print(f"{'i':3} {'V':>10} {'Vn':>10} {'Vp':>10} {'W (nm)':>10} {'N(W) (cm-3)':>10} {'N_meas(W) (cm-3)':>10} {'Wn (nm)':>10} {'Wp (nm)':>10} {'C (F/m2)':>10} {'1/C^2 (m^2F^-2)':>10} {'Wmeas (nm)':>10} {'Cmeas (F/m2)':>10} {'1/Cmeas^2 (m4F^-2)':>10}")
    else:
        print(f"{'i':3} {'V':>10} {'Vn':>10} {'Vp':>10} {'W (nm)':>10} {'N(W) (cm-3)':>10} {'Wn (nm)':>10} {'Wp (nm)':>10} {'C (F/m2)':>10} {'1/C^2 (m^2F^-2)':>10}")
    for i in range(nV):
        _W, _N = pnj.cal_N_W(i, V, C2inv, Wn, Wp)
        N.append(_N)
        if C_meas:
            _W, _N = pnj.cal_N_W(i, V, Cm2inv_meas, Wn, Wp)
            N_meas.append(_N)
            W_meas.append(_W)

        if i % nskip == 0:
            if C_meas:
                print(f"{i:3} {V[i]:10.4g} {Vn[i]:10.4g} {Vp[i]:10.4g} {W[i]*1e9:10.4g} {N[i]*1.0e-6:10.4g} {N_meas[i]*1.0e-6:10.4g} {Wn[i]*1e9:10.4g} {Wp[i]*1e9:10.4g} {C[i]:10.4g} {C2inv[i]:10.4g} {_W*1e9:10.4g} {Cm2_meas[i]:10.4g} {Cm2inv_meas[i]:10.4g}")
            else:
                print(f"{i:3} {V[i]:10.4g} {Vn[i]:10.4g} {Vp[i]:10.4g} {W[i]*1e9:10.4g} {N[i]*1.0e-6:10.4g} {Wn[i]*1e9:10.4g} {Wp[i]*1e9:10.4g} {C[i]:10.4g} {C2inv[i]:10.4g}")
            
    print()
    print(f"V at full depletion in n layer: Vn_fdep={Vn_fdep:.4g} V")

    S = cparams.S
    print()
    print(f"Save data to {cparams.outfile}")
    if C_meas:
        df = pd.DataFrame({'V (V)': V, 'W (nm)': np.array(W) * 1.0e9, 'N (cm^-3)': np.array(N) * 1.0e-6, 'N_meas (cm^-3)': np.array(N_meas) * 1.0e-6, 
                       'Wn (m)': np.array(Wn) * 1.0e9, 'Wp (m)': np.array(Wp) * 1.0e9,
                       'C (F/m^2)': C, '1/C2 (m^4/F^-2)': C2inv, 
                       'C (F)': np.array(C)*S, '1/C2 (F^-2)': np.array(C2inv)/S/S, 
                       'Cmeas (F/m2)': Cm2_meas, '1/Cmeas2 (m^4F^-2)': Cm2inv_meas})
    else:
        df = pd.DataFrame({'V (V)': V, 'W (nm)': np.array(W) * 1.0e9, 'N (cm^-3)': np.array(N) * 1.0e-6, 
                       'Wn (m)': np.array(Wn) * 1.0e9, 'Wp (m)': np.array(Wp) * 1.0e9,
                       'C (F/m^2)': C, '1/C2 (m^4/F^-2)': C2inv, 
                       'C (F)': np.array(C)*S, '1/C2 (F^-2)': np.array(C2inv)/S/S})
    df.to_excel(cparams.outfile, index = False)

#=========================================
# plot
#=========================================
    if not cparams.plot:
        return
        
    plt.rcParams['font.size'] = cparams.fontsize
    fig = plt.figure(figsize = (12, 8))
    ax1 = fig.add_subplot(2, 3, 1)
    ax2 = fig.add_subplot(2, 3, 2)
    ax3 = fig.add_subplot(2, 3, 3)
    ax4 = fig.add_subplot(2, 3, 4)
    ax5 = fig.add_subplot(2, 3, 5)  
    ax6 = fig.add_subplot(2, 3, 6)  

    ax1.plot(V, np.array(C2inv),       label = '1/$C^2$', linestyle = '', marker = 'o', markersize = 5.0)
    if C_meas:
        ax1.plot(V, np.array(Cm2inv_meas), label = '1/$C_{meas}^2$', color = 'black')
    ax1.set_xlabel('V (V)', fontsize = cparams.fontsize)
    ax1.set_ylabel('1/$C^2$ (m$^4$F$^{-2}$)',        fontsize = cparams.fontsize)
#    ax1.set_ylabel('1/$C_{meas}^2$ (m$^4$F$^{-2}$)', fontsize = cparams.fontsize)
#    ax1.set_yscale('log')
    ax1.legend(fontsize = cparams.legend_fontsize)

    ax2.plot(V, np.array(C),        label = 'C', linestyle = '', marker = 'o', markersize = 5.0)
    if C_meas:
        ax2.plot(V, np.array(Cm2_meas), label = 'C$_{meas}$', color = 'black')
    ax2.set_xlabel('V (V)', fontsize = cparams.fontsize)
    ax2.set_ylabel('C (F)', fontsize = cparams.fontsize)
#    ax2.set_yscale('log')
    ax2.legend(fontsize = cparams.legend_fontsize)

    colors = []
    for i in range(len(Wn)):
        if Wn[i] >= cparams.dn and Wp[i] >= cparams.dp:
            colors.append('black')
        elif Wn[i] >= cparams.dn:
            colors.append('red')
        else:
            colors.append('blue')

    ax3.plot(V, np.array(Wn) * 1.0e9, label = '$W_n$ (nm)', color = 'red', linewidth = 0.5)
    ax3.scatter(V, np.array(Wn) * 1.0e9, color = colors, marker = '>', s = 3.0)
    ax3.plot(V, np.array(Wp) * 1.0e9, label = '$W_p$ (nm)', color = 'blue', linewidth = 0.5)
    ax3.scatter(V, np.array(Wp) * 1.0e9, color = colors, marker = '^', s = 3.0)
    ax3.plot(V, np.array(W) * 1.0e9,  label = '$W$ (nm)',   color = 'black', linewidth = 1.0)
    ax3.scatter(V, np.array(W) * 1.0e9, color = colors, marker = 'o', s = 6.0)
    ax3.set_xlabel('V (V)', fontsize = cparams.fontsize)
    ax3.set_ylabel('W (nm)', fontsize = cparams.fontsize)
    ax3.set_yscale('log')
    ax3.legend(fontsize = cparams.legend_fontsize)

    if C_meas:
        wnm = np.array(W_meas) * 1.0e9
        ax4.plot(wnm, np.array(N) * 1.0e-6, label = '$N$ (cm$^{-3}$)')
        ax4.plot(wnm, np.array(N_meas) * 1.0e-6, label = '$N_{meas}$ (cm$^{-3}$)')
        ax4.set_xlabel('W (nm)',        fontsize = cparams.fontsize)
    else:
        ax4.plot(V, np.array(N) * 1.0e-6, label = '$N$ (cm$^{-3}$)')
        ax4.plot(V, np.array(N_meas) * 1.0e-6, label = '$N_{meas}$ (cm$^{-3}$)')
        ax4.set_xlabel('V (V)',        fontsize = cparams.fontsize)
    ax4.set_ylabel('N (cm$^{-3}$)', fontsize = cparams.fontsize)
    ax4.set_yscale('log')
    ax4.legend(fontsize = cparams.legend_fontsize)

    V_Cp_fdep = []
    V_Cp_ndep = []
    Cp2inv_fdep = []
    Cp2inv_ndep = []
    for iV, _V in enumerate(V_Cp):
        if _V < Vn_fdep:
            V_Cp_fdep.append(_V)
            Cp2inv_fdep.append(Cp2inv[iV])
        else:
            V_Cp_ndep.append(_V)
            Cp2inv_ndep.append(Cp2inv[iV])

    ax5.plot(V_Cp_fdep, np.array(Cp2inv_fdep), label = 'n full depleted 1.0 / $C_p^2$ (m$^4$F$^{-2}$)', color = 'red', linewidth = 0.5)
    ax5.plot(V_Cp_ndep, np.array(Cp2inv_ndep), label = 'n not depleted 1.0 / $C_p^2$ (m$^4$F$^{-2}$)', color = 'black', linewidth = 0.5)
    ax5.set_xlabel('V (V)', fontsize = cparams.fontsize)
    ax5.set_ylabel('1/$C_p^2$ (m$^4$F$^{-2}$)',        fontsize = cparams.fontsize)

    ax6.plot(np.array(W_Cp) * 1.0e9, np.array(N_Cp) * 1.0e-6, label = '$N$ in p (cm$^{-3}$)')
    ax6.set_xlabel('W (nm)',        fontsize = cparams.fontsize)
    ax6.set_ylabel('N (cm$^{-3}$)', fontsize = cparams.fontsize)
    ax6.set_yscale('log')
    ax6.legend(fontsize = cparams.legend_fontsize)


    plt.tight_layout()

    plt.pause(0.0001)
    input("\nPress ENTER to terminate>>")


def main():
    cparams = initialize()
    update_variables(cparams)
    
    if cparams.mode == 'NW':
        plot_NW(cparams)
    else:
        print(f"\nError in main(): Invalide mode [{cparams.mode}]\n")
        input("Press ENTER to terminate>>")
        exit()

    exit()


if __name__ == "__main__":
    main()