#!/usr/bin env python3
from formulary import *
from scipy.constants import Boltzmann as kB, hbar, electron_volt


# DEVICES
# metal sc: schottky barrier
def schottky_barrier():
    fig, axs = plt.subplots(1, 3, figsize=(width_formula, height_default*0.6))
    WD = 5
    q = 1
    ND = 1
    eps = 11
    dx = WD/10
    xs = np.linspace(-WD/5, 6/5*WD, 300)
    rho_S = q*ND
    Q = rho_S * WD
    rho_M = -Q/dx
    @np.vectorize
    def rho_approx(x):
        if x < -dx: return 0.0
        if x < 0: return rho_M
        if x < WD: return rho_S
        return 0.0
    rhos_approx = rho_approx(xs)

    @np.vectorize
    def E(x):
        if x < -dx: return 0.0
        if x < 0: return -rho_M/eps * (-dx-x)
        if x < WD: return -rho_S/eps * (WD-x)
        return 0.0
    Es = E(xs)

    @np.vectorize
    def phi(x):
        # if x < -dx: return 0.0
        # if x < 0: return -rho_M/(2*eps) * (dx**2-(dx-x)**2)
        if x < 0: return 0.0
        if x < WD: return rho_S/(2*eps) * (WD**2-(WD-x)**2)
        return q*ND/(2*eps) * WD**2
    phis = phi(xs)

    for ax in axs:
        ax.set_xlabel("$x$")
        ax.set_xticks([0,WD])
        ax.set_xticklabels(["0", r"$W_\text{D}$"])
        ax.set_yticks([0])
        ax.set_yticklabels(["0"])
    axs[0].plot(xs, rhos_approx)
    axs[0].set_ylabel(r"$\rho(x)$")
    axs[1].plot(xs, Es)
    axs[1].set_ylabel(r"$\mathcal{E}(x)$")

    axs[2].plot(xs, phis)
    axs[2].set_ylabel(r"$\phi(x)$")

    return fig





# Charge carrier density
def fn_i(T, NV, NC, Egap):
    return np.sqrt(NV*NC) * np.exp(-Egap*electron_volt/(2*kB*T))

def fn_freeze(T, NV, NC, Egap):
    return np.sqrt(NV*NC) * np.exp(-0.07*Egap*electron_volt/(2*kB*T))

def charge_carrier_density():
    fig, ax = plt.subplots(1, 1, figsize=size_formula_normal_default)
    # Ts = np.power(10, np.linspace(-4, 0, 100))
    N = 500
    ts = np.linspace(0.01, 20, N)
    Ts = 1000/ts
    # Ts = np.linspace(2000, 50, N)
    # ts = 1000/Ts
    # ts = 1/Ts

    NV = 1e23
    NC = 1e21
    Egap = 2.4

    n_is = fn_i(Ts, NV, NC, Egap)

    n_sat = np.empty(Ts.shape)
    n_sat.fill(1e15)


    idx1 = np.argmin(np.abs(n_is-n_sat))
    print(idx1, N)
    # TODO make quadratic simple
    n_total = blend_arrays(n_is, n_sat, idx1, idx1+10, "linear")

    n_freeze = fn_freeze(Ts, NV, NC, Egap)
    idx2 = np.argmin(np.abs(n_freeze-n_sat))

    print(idx1, idx2, N)
    n_total = blend_arrays(n_total, n_freeze, idx2-10, idx2+10, "quadratic_simple")

    # ax.plot(ts, n_is, label="Intrinsic")
    # ax.plot(ts, n_sat, label="Saturation")
    # ax.plot(ts, n_freeze, label="Freeze-out")
    # ax.plot(ts, n_total, label="Total", linestyle="dashed", color="black")
    n_total2 = n_is.copy()
    n_total2[idx1:idx2] = n_sat[idx1:idx2]
    n_total2[idx2:] = n_freeze[idx2:]
    ax.plot(ts, n_is, label="Intrinsic", linestyle="dashed")
    ax.plot(ts, n_total2, label="Total", color="black")
    ax.set_yscale("log")
    ax.set_ylim(1e13, 1e17)

    idx = int(N*0.9)
    ax.text(ts[idx], n_freeze[idx], "Freeze-out", ha="right", va="top")
    idx = int(N*0.5)
    ax.text(ts[idx], n_sat[idx], "Saturation", ha="center", va="bottom")
    idx = np.argmin(np.abs(n_is-3e16))
    ax.text(ts[idx+10], n_is[idx], "Intrinsic", ha="left", va="bottom")
    # ax.set_xlim(ts[0], ts[idx1+N//6])
    # ax.legend()
    ax.set_xlabel(r"$1000/T\,[\si{\per\kelvin}]$")
    ax.set_ylabel(r"$n\,[\si{\per\cm^3}]$")


    return fig


def test():
    fig, ax = plt.subplots()
    N = 100
    left = np.empty(N)
    left.fill(4.0)
    left = np.linspace(5.0, 0.0, N)
    right = np.empty(N)
    right.fill(2.5)
    right = np.linspace(3.0, 2.0, N)
    ax.plot(left, label="l",linestyle="dashed")
    ax.plot(right, label="r",linestyle="dashed")

    total_lin = blend_arrays(left, right, 40, 60, "linear")
    ax.plot(total_lin, label="lin")
    # total_tanh = blend_arrays(left, right, 40, 60, "tanh")
    # ax.plot(total_tanh)
    total_q = blend_arrays(left, right, 40, 60, "quadratic_simple")
    ax.plot(total_q, label="q")
    ax.legend()
    return fig






if __name__ == '__main__':
    # export(test(), "cm_sc_charge_carrier_density")
    export(schottky_barrier(), "cm_sc_devices_metal-n-sc_schottky")
    # export(charge_carrier_density(), "cm_sc_charge_carrier_density")