#!/usr/bin env python3
from formulary import *
from scipy.constants import gas_constant, Avogadro, elementary_charge

Faraday = Avogadro * elementary_charge

@np.vectorize
def fbutler_volmer_anode(ac, z, eta, T):
    return np.exp((1-ac)*z*Faraday*eta/(gas_constant*T))

@np.vectorize
def fbutler_volmer_cathode(ac, z, eta, T):
    return -np.exp(-ac*z*Faraday*eta/(gas_constant*T))

def fbutler_volmer(ac, z, eta, T):
    return fbutler_volmer_anode(ac, z, eta, T) + fbutler_volmer_cathode(ac, z, eta, T)

def butler_volmer():
    fig, ax = plt.subplots(figsize=size_half_third)
    ax.set_xlabel("$\\eta$ [V]")
    ax.set_ylabel("$j/j_0$")
    etas = np.linspace(-0.1, 0.1, 400)
    T = 300
    z = 1.0
    # other a
    ac2, ac3 = 0.2, 0.8
    i2 = fbutler_volmer(0.2, z, etas, T)
    i3 = fbutler_volmer(0.8, z, etas, T)
    ax.plot(etas, i2, color="blue", linestyle="dashed", label=f"$\\alpha_\\text{{C}}={ac2}$")
    ax.plot(etas, i3, color="green", linestyle="dashed", label=f"$\\alpha_\\text{{C}}={ac3}$")
    # 0.5
    ac = 0.5
    irel_anode = fbutler_volmer_anode(ac, z, etas, T)
    irel_cathode = fbutler_volmer_cathode(ac, z, etas, T)
    ax.plot(etas, irel_anode,  color="gray")
    ax.plot(etas, irel_cathode, color="gray")
    ax.plot(etas, irel_cathode + irel_anode, color="black", label=f"$\\alpha_\\text{{C}}=0.5$")
    ax.grid()
    ax.legend()
    ylim = 6
    ax.set_ylim(-ylim, ylim)
    return fig

@np.vectorize
def ftafel_anode(ac, z, eta, T):
    return 10**((1-ac)*z*Faraday*eta/(gas_constant*T*np.log(10)))

@np.vectorize
def ftafel_cathode(ac, z, eta, T):
    return -10**(-ac*z*Faraday*eta/(gas_constant*T*np.log(10)))

def tafel():
    i0 = 1
    ac = 0.2
    z = 1
    T = 300
    eta_max = 0.2
    etas = np.linspace(-eta_max, eta_max, 400)
    i = np.abs(fbutler_volmer(ac, z, etas ,T))
    iright = i0 * np.abs(ftafel_cathode(ac, z, etas, T))
    ileft = i0 * ftafel_anode(ac, z, etas, T)

    fig, ax = plt.subplots(figsize=size_half_third)
    ax.set_xlabel("$\\eta$ [V]")
    ax.set_ylabel("$\\log_{10}\\left(\\frac{|j|}{j_0}\\right)$")
    # ax.set_ylabel("$\\log_{10}\\left(|j|/j_0\\right)$")
    ax.set_yscale("log")
    # ax.plot(etas, linear, label="Tafel slope")
    ax.plot(etas[etas >= 0], ileft[etas >= 0],  linestyle="dashed", color="gray", label="Tafel Approximation")
    ax.plot(etas[etas <= 0], iright[etas <= 0], linestyle="dashed", color="gray")
    ax.plot(etas, i, label=f"Butler-Volmer $\\alpha_\\text{{C}}={ac:.1f}$")
    ax.legend()
    ax.grid()
    return fig

if __name__ == '__main__':
    export(butler_volmer(), "ch_butler_volmer")
    export(tafel(), "ch_tafel")