import numpy as np
import matplotlib.pyplot as plt

a = 0.5
F = 96485
R = 8.314
T = 298

CR = 1
CO = 1

k0 = 1e-5

E_rd = np.linspace(-500e-3, -400e-3, 10)
E_ox = np.linspace(400e-3, 500e-3, 10)

i_ox = lambda E: F * k0 * CR * np.exp(((1-a)*F*E)/(R*T))
i_rd = lambda E: -F * k0 * CO * np.exp(((-a)*F*E)/(R*T))

i = lambda E: i_ox(E) + i_rd(E)

tafel_i = lambda E: np.log(np.abs(i(E)))

koeff_ox = np.polyfit(E_ox, tafel_i(E_ox), 1)
koeff_rd = np.polyfit(E_rd, tafel_i(E_rd), 1)

E_akse = np.linspace(-500e-3, 500e-3, 100)

fig, axs = plt.subplots(1,2, figsize = [10,5])

plt.sca(axs[0])
plt.plot(E_rd, i_rd(E_rd), color = 'red')
plt.plot(E_ox, i_ox(E_ox), color = 'blue')
plt.plot(E_akse, i(E_akse), color = 'green', linestyle = '--')
plt.xlabel(r'$E-E_n$')
plt.ylabel('i')

plt.sca(axs[1])
plt.plot(E_akse, np.log(np.abs(i(E_akse))), color = 'green', linestyle = '--')
plt.plot(E_akse[:55], np.polyval(koeff_rd, E_akse[:55]), color = 'red')
plt.plot(E_akse[45:], np.polyval(koeff_ox, E_akse[45:]), color = 'blue')

plt.plot([E_akse[0], 0], [koeff_rd[1], koeff_rd[1]], color = 'black',
         linestyle = '--', label = r'$\ln{i_n} = $'+str(round(koeff_rd[1],4)))
plt.legend()
plt.ylabel(r'$\ln{| i |}$')
plt.xlabel(r'$E-E_n$')

plt.savefig('Finding_exchange_current', dpi = 600)

plt.show()