# import external modules
import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.size'] = 18
import scipy.signal as signal
import os

# import Comphy modules
from ComphyV3.CoreV3.MOS_1D import MOS_1D
from ComphyV3.CoreV3.Physics import BandGapModel_Tdep1
from ComphyV3.CoreV3.Physics import CarrierModel_Electrons_1
from ComphyV3.CoreV3.Physics import CarrierModel_Holes_1
from ComphyV3.CoreV3.Channel import Channel
from ComphyV3.CoreV3.InsulatingLayer import InsulatingLayer
from ComphyV3.CoreV3.Metal import Metal
from ComphyV3.CoreV3.TrapBand import TrapBand_2S_NMP
from ComphyV3.CoreV3.TrapBand import GeneratedInterfaceTraps
from ComphyV3.CoreV3.GridSampler import GaussianGridSampler
from ComphyV3.CoreV3.GridSampler import InterfaceGridSampler

from ComphyV3.Utils.Plotter import plot_band_diagram, plot_potential_profile, plot_electric_field, plot_surface_potential


# create new output directory
output_dir = './Results/'
os.makedirs(output_dir, exist_ok=True)


########################################################################################################################
# Basic Usage Part1: Create your own device
########################################################################################################################

Temperature = 300.0
my_device = MOS_1D(
    Tinit=Temperature,  # initial temperature in K
)

######################################################################
# Step 2: Add semiconductor channel:

bandgap_model = BandGapModel_Tdep1(
    Eg0=1.206,  # band gap parameter in eV
    Eg1=-2.73E-4,  # band gap parameter in eV/K
)

e_model = CarrierModel_Electrons_1(
    Ncv0=2.541E25,  # effective density of states of conduction band
    Mc=6,  # equivalent minima of conduction band
    mt0=0.1905,  # transversal effective mass
    ml=0.9163,  # longitudinal effective mass
    bgmodel=bandgap_model,  # bandgap model
)

h_model = CarrierModel_Holes_1(
    Ncv0=2.54E25,  # effective density of states of valence band
    mEff_a=0.4435870,  # effective mass parameter
    mEff_b=0.3609528e-2,  # effective mass parameter in K⁻¹
    mEff_c=0.1173515e-3,  # effective mass parameter in K⁻²
    mEff_d=0.1263218e-5,  # effective mass parameter in K⁻³
    mEff_e=0.3025581e-8,  # effective mass parameter in K⁻⁴
    mEff_f=0.4683382e-2,  # effective mass parameter in K⁻¹
    mEff_g=0.2286895e-3,  # effective mass parameter in K⁻²
    mEff_h=0.7469271e-6,  # effective mass parameter in K⁻⁶
    mEff_i=0.1727481e-8,  # effective mass parameter in K⁻⁸
)

Si = Channel(
    Nd=1.4E23,  # donor concentration in m⁻³
    Na=1.0E16,  # acceptor concentration mass in m⁻³
    name='Si',  # material name
    Eg=bandgap_model,  # band gap model of channel
    m_neff=e_model,  # electron effective mass (or carrier model)
    m_peff=h_model,  # hole effective mass (or carrier model)
    Xsi=4.05,  # electron affinity in eV
    kval=11.68,  # relative permittivity
    ccs_e=2.0E-23,  # capture cross section for electrons in m⁻²
    ccs_h=2.0E-23,  # capture cross section for holes in m⁻²
)
my_device.channel = Si

######################################################################
# Step 3: Add oxide layers:

SiO2 = InsulatingLayer(
    thickness=1.0E-9,  # thickness of layer in m
    kval=3.9,  # relative permittivity
    Eg=9.0,  # band gap in eV (or bandgap model)
    Xsi=0.82,  # electron affinity in eV
    material_name='SiO2',  # name of layer-material
    mt=0.35,  # relative tunneling mass
)
my_device.insulating_layers.append(SiO2)

HfO2 = InsulatingLayer(
    thickness=2.0E-9,  # thickness of layer in m
    kval=20.0,  # relative permittivity
    Eg=5.8,  # band gap in eV (or bandgap model)
    Xsi=1.972,  # electron affinity in eV
    material_name='HfO2',  # name of layer-material
    mt=0.17,  # relative tunneling mass
)
my_device.insulating_layers.append(HfO2)


######################################################################
# Step 4: Add gate:

WF = 4.05 + 0.56 + 0.1
gate = Metal(
    workfunction=WF,  # work function in eV
    name='metal',  # name of gate-material
    ccs=2.0E-23,  # capture cross section for electrons in m⁻²
)
my_device.gate = gate

######################################################################
# Step 5: Print basic device properties:

# print basic device properties

print("Vfb0 = {:.3f} V".format(my_device.get_Vfb0(T=Temperature)))

print("Vth0 = {:.3f} V".format(my_device.get_Vth0(T=Temperature)))

print("EOT = {:.3f} nm".format(my_device.EOT * 1.0E9))

print("Nd = {:.3e} cm⁻³".format(my_device.channel.Nd * 1.0E-6))

print("Na = {:.3e} cm⁻³".format(my_device.channel.Nd * 1.0E-6))

print("Ef = {:.3f} eV".format(my_device.channel.get_Fermilevel(T=Temperature)))

for i, layer in enumerate(my_device.insulating_layers):
    print("layer {}:".format(i), layer)

######################################################################
# Step 5: Plot flatband diagram:

# initialize devices
my_device.initialize(Vg=0.521, T=Temperature)

# plot band diagram
plot_band_diagram(my_device, file=output_dir + "flat_band_diagram.png")


########################################################################################################################
# Basic Usage Part2: Adding charge traps to your device
########################################################################################################################


######################################################################
# Step 1: Modeling the recoverable BTI component

grid_sampler = GaussianGridSampler(
    dx=1.0E-10,  # grid spacing for sampling xt in m
    dEt=0.05,  # grid spacing for sampling of Et in eV
    dEr=0.1,  # grid spacing for sampling of S in eV
    non_coverage=1.0E-4,  # sets boundary of gaussian grid
)

SiO2_shallow_trap_band = TrapBand_2S_NMP(
    x_start=0.0,  # lower bound of spatial defect distribution in m
    x_end=SiO2.thickness,  # upper bound of spatial defect distribution in m
    Et=-3.48,  # mean energy level in eV
    Et_sigma=0.15,  # standard dev. of energy level in eV
    Er=3.82,  # mean relaxation energy in eV
    Er_sigma=1.36,  # standard dev. of relaxation energy in eV
    R=0.407,  # curvature ratio
    Nt=1.5e25,  # defect concentration in m⁻³
    trap_type='acceptor',  # trap type
    sampler=grid_sampler,  # a appropriate grid sampler
    device=my_device,  # the device to which you want to add the traps
)

my_device.trap_bands['SiO2_shallow_trap_band'] = SiO2_shallow_trap_band

SiO2_deep_trap_band = TrapBand_2S_NMP(
    x_start=0.0,  # lower bound of spatial defect distribution in m
    x_end=SiO2.thickness,  # upper bound of spatial defect distribution in m
    Et=-5.65,  # mean energy level in eV
    Et_sigma=0.20515,  # standard dev. of energy level in eV
    Er=9.3438,  # mean relaxation energy in eV
    Er_sigma=2.5526,  # standard dev. of relaxation energy in eV
    R=1.809,  # curvature ratio
    Nt=4.449E26,  # defect concentration in m⁻³
    trap_type='donor',  # trap type
    sampler=grid_sampler,  # a appropriate grid sampler
    device=my_device,  # the device to which you want to add the traps
)

my_device.trap_bands['SiO2_deep_trap_band'] = SiO2_deep_trap_band

HfO2_shallow_trap_band = TrapBand_2S_NMP(
    x_start=1.0E-9,  # lower bound of spatial defect distribution in m
    x_end=3.0E-9,  # upper bound of spatial defect distribution in m
    Et=-3.18,  # mean energy level in eV
    Et_sigma=0.25455,  # standard dev. of energy level in eV
    Er=0.8684,  # mean relaxation energy in eV
    Er_sigma=0.4045,  # standard dev. of relaxation energy in eV
    R=0.0018121,  # curvature ratio
    Nt=9.0534E25,  # defect concentration in m⁻³
    trap_type='acceptor',  # trap type
    sampler=grid_sampler,  # a appropriate grid sampler
    device=my_device,  # the device to which you want to add the traps
)

my_device.trap_bands['HfO2_shallow_trap_band'] = HfO2_shallow_trap_band

HfO2_deep_trap_band = TrapBand_2S_NMP(
    x_start=1.0E-9,  # lower bound of spatial defect distribution in m
    x_end=3.0E-9,  # upper bound of spatial defect distribution in m
    Et=-5.232,  # mean energy level in eV
    Et_sigma=0.23203,  # standard dev. of energy level in eV
    Er=9.2285,  # mean relaxation energy in eV
    Er_sigma=1.4948,  # standard dev. of relaxation energy in eV
    R=1.5287,  # curvature ratio
    Nt=4.1597E26,  # defect concentration in m⁻³
    trap_type='donor',  # trap type
    sampler=grid_sampler,  # a appropriate grid sampler
    device=my_device,  # the device to which you want to add the traps

)

my_device.trap_bands['HfO2_deep_trap_band'] = HfO2_deep_trap_band


######################################################################
# Step 2: Modeling the quasi-permanent BTI component

interface_grid_sampler = InterfaceGridSampler(
    d_eps1=0.1,  # grid spacing for sampling eps_1 in eV
    d_eps2=0.001,  # grid spacing for sampling eps_2 in eV
)

interface_traps = GeneratedInterfaceTraps(
    eps1=2.5297,  # mean activation energy barrier in eV
    eps1_sigma=0.92612,  # standard dev. of activation energy barriers in eV
    eps2=1.649,  # mean passivation energy barrier in eV
    eps2_sigma=0.060114,  # standard dev. of passivation energy barriers in eV
    k0=1.2336E13,  # transition rate for zero barrier in s⁻¹
    gamma=1.7052e-10,  # field dependence of activation barrier in eV m V⁻¹
    Nt=1.9096E17,  # curvature ratio
    trap_type='donor',  # trap type
    sampler=interface_grid_sampler,  # interface grid sampler
    device=my_device,  # the device to which you want to add the traps
)



######################################################################
# Step 3: Initializing the device
my_device.initialize(
    Vg=-0.5,  # initial gate voltage in V
    T=300.0,  # initial temperature in K
    record_band_occus=False,  # indicates whether or not to save all occupancies of all bands at all timesteps
    debug=True, threshold=1.0E-65
)


plot_band_diagram(my_device, file=output_dir + "band_diagram_equilibrium.png", draw_samples=None)


########################################################################################################################
# Basic Usage Part3: Performing a simple BTI simulation
########################################################################################################################

######################################################################
# Step 1: Print basic device properties

print("Vfb0 = {:.3f} V".format(my_device.get_Vfb0(T=Temperature)))

print("Vth0 = {:.3f} V".format(my_device.get_Vth0(T=Temperature)))

print("Vfb = {:.3f} V".format(my_device.Vfb))

print("Vth = {:.3f} V".format(my_device.Vth))

######################################################################
# Step 2: Plot semiconductor potential

plot_surface_potential(my_device, file=output_dir + "surface_potential.png",)

######################################################################
# Step 3: Calculate potential profile and electric field across the stack

plot_potential_profile(my_device, file=output_dir + "potential_profile.png")

plot_electric_field(my_device, file=output_dir + "electric_field.png")

######################################################################
# Step 4: Define stress pattern

# define time steps
time_per_step = 1
number_steps = 250
time = np.linspace(0.0, time_per_step * number_steps, number_steps)

# define constant temperature profile
temperature = np.full_like(time, 300.0)

# define gate voltage profile
min = -2.0
max = -0.5
frequency = 1.0 / 40.0
Vg = min + (max - min)*(signal.square(2 * np.pi * frequency * time)/2.0 + 0.5)

# plot gate voltage profile
figure = plt.figure(figsize=(8.0, 6.0))
plt.plot(time, Vg)
plt.plot(time, Vg, 'o')
plt.title("stress waveform")
plt.xlabel("time [s]")
plt.ylabel("gate voltage [V]")
plt.tight_layout()
plt.grid()
plt.savefig(output_dir + "stress_waveform.png")
plt.show()


######################################################################
# Step 5: Calculate shift of threshold voltage

# apply stress pattern to device

# let's initialize the device with the initial voltage and temperature
my_device.initialize(Vg=Vg[0], T=temperature[0])

# Now, let's retrieve the threshold voltage of the unstressed device as reference
Vth_ref = my_device.Vth

# Now, we create a list to collect the threshold voltage shifts
threshold_voltage_shifts = [0.0]

# Now, we loop over the scenario we have set out.
# We skip the first point as that is our initialization point.
for idx in range(1, len(time)):

    # here we actually apply the voltage for a given timestep
    my_device.apply_gate_voltage(
        Vg=Vg[idx],
        T=temperature[idx],
        dt=time[idx] - time[idx-1],
    )
    # now let's retrieve the threshold voltage after application of stress
    new_Vth = my_device.Vth

    # finally, we store the threshold voltage w.r.t to the reference point
    threshold_voltage_shifts.append(new_Vth - Vth_ref)


# Now that we have the data, let's plot it
figure = plt.figure(figsize=(8.0, 6.0))
plt.plot(time, np.abs(threshold_voltage_shifts))
plt.title("BTI sequence")
plt.xlabel("time [s]")
plt.ylabel("threshold voltage shift [V]")
plt.tight_layout()
plt.grid()
plt.savefig(output_dir + "bti_sequence.png")
plt.show()
