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

Faraday = Avogadro * elementary_charge

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

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

def fbutler_volmer(ac, z, eta, T):
    return fbutler_volmer_left(ac, z, eta, T) + fbutler_volmer_right(ac, z, eta, T)

def butler_volmer():
    fig, ax = plt.subplots(figsize=size_half_third)
    ax.set_xlabel("$\\eta$")
    ax.set_ylabel("$i/i_0$")
    etas = np.linspace(-0.1, 0.1, 400)
    T = 300
    z = 1.0
    # other a
    alpha2, alpha3 = 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={alpha2}$")
    ax.plot(etas, i3, color="green", linestyle="dashed", label=f"$\\alpha={alpha3}$")
    # 0.5
    ac = 0.5
    irel_left = fbutler_volmer_left(ac, z, etas, T)
    irel_right = fbutler_volmer_right(ac, z, etas, T)
    ax.plot(etas, irel_left,  color="gray")
    ax.plot(etas, irel_right, color="gray")
    ax.plot(etas, irel_right + irel_left, color="black", label=f"$\\alpha=0.5$")
    ax.grid()
    ax.legend()
    ylim = 6
    ax.set_ylim(-ylim, ylim)
    return fig

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