Tafel plots (Atkins 8. Ed. Ch. 25), using twin x axis with twiny

WBzT_current_voltage.py

import matplotlib.pyplot as plt
from numpy import array, ones, log, exp, concatenate, linspace, mean, sqrt
from numpy.linalg import inv

eta=array([float(x) for x in '50 100 150 200 250'.split(' ')]) # mV
i=array([float(x) for x in '8.8 25.0 58.0 131 298'.split(' ')]) # mA
j=i/2 # mA/cm^2

ln_j=log(j)

ln_j_fit=ln_j[1:]
eta_fit=eta[1:]

h=array([ones(len(eta_fit)),eta_fit]).T # designmatrix

pseudoinverse=inv((h.T.dot(h))).dot(h.T)
params=pseudoinverse.dot(ln_j_fit)

h1=concatenate([[[1,50]],h])
alpha=1-params[1]*1000*8.3145*298.15/96485

plt.subplot(2,1,1)
plt.plot(eta,ln_j,'o',eta,h1.dot(params))
plt.ylabel(r'$ln\left(j/(mA cm^{-2})\right)$')
plt.xlabel(r'$\eta / mV$')
plt.text(150,2,
    '$j_0='+'{:0.5}'.format(exp(params[0]))+r'~mA\cdot cm^{-2}'+'$'+'\n'
    r'$\frac{(1-\alpha)F n}{R T}='+
    '{:0.5}'.format(params[1]*1000)+r'\frac{1}{V}\Rightarrow\alpha='+'{:0.5}'.format(alpha)+
    '$',ha='left')
plt.title('Atkins 8. E. Ex. 25.4')

plt.subplot(2,1,2)
eta_compl=linspace(-300,300,40)/1000 # V
j0=exp(params[0])
j_compl=j0*(exp((1-alpha)*96485*1*eta_compl/(8.3145*298.15))-exp(-alpha*96485*1*eta_compl/(8.3145*298.15)))
low_eta=array([-0.12,0.12])
j_low_eta=j0*low_eta*95485*1/(8.3145*298.15)
plt.plot(eta/1000,j,'o',eta_compl,j_compl,'-',low_eta,j_low_eta,'--')
plt.xlabel(r'$\eta / V$')
plt.ylabel('$j/(mA cm^{-2})$')
plt.legend(['exp','Butler-Volmer',r'$-0,12~V\leq\eta\leq 0,12~V$'])
plt.tight_layout()

eta=array([float(x) for x in '−50 −100 −150 −200 −250 −300'.replace('−','-').split(' ')]) # mV
i=array([float(x) for x in '−0.3 −1.5 −6.4 −27.6 −118.6 −510'.replace('−','-').split(' ')]) # mA
j=i/2 # mA/cm^2

ln_minus_j=log(-j)

ln_minus_j_fit=ln_minus_j[1:]
eta_fit=eta[1:]

h=array([ones(len(eta_fit)),eta_fit]).T # designmatrix

pseudoinverse=inv((h.T.dot(h))).dot(h.T)
params=pseudoinverse.dot(ln_minus_j_fit)

h1=concatenate([[[1,-50]],h])
alpha=-params[1]*1000*8.3145*298.15/96485

plt.figure()
plt.subplot(2,1,1)
plt.plot(eta,ln_minus_j,'o',eta,h1.dot(params))
plt.ylabel(r'$ln\left(j/(mA cm^{-2})\right)$')
plt.xlabel(r'$\eta / mV$')
plt.text(-300,0,
    '$j_0='+'{:0.5}'.format(exp(params[0]))+r'~mA\cdot cm^{-2}'+'$'+'\n'
    r'$\frac{-\alpha F n}{R T}='+
    '{:0.5}'.format(params[1]*1000)+r'\frac{1}{V}\Rightarrow\alpha='+'{:0.5}'.format(alpha)+
    '$',ha='left')
plt.title('Atkins 8. E. Self test 25.7')

plt.subplot(2,1,2)
eta_compl=linspace(-300,700,40)/1000 # V
j0=exp(params[0])
j_compl=j0*(exp((1-alpha)*96485*1*eta_compl/(8.3145*298.15))-exp(-alpha*96485*1*eta_compl/(8.3145*298.15)))
low_eta=array([-0.12,0.12])
j_low_eta=j0*low_eta*95485*1/(8.3145*298.15)
plt.plot(eta/1000,j,'o',eta_compl,j_compl,'-',low_eta,j_low_eta,'--')
plt.xlabel(r'$\eta / V$')
plt.ylabel('$j/(mA cm^{-2})$')
plt.legend(['exp','Butler-Volmer',r'$-0,12~V\leq\eta\leq 0,12~V$'])
plt.tight_layout()


e=array([float(x) for x in '0.388 0.365 0.350 0.335'.split(' ')])*-1*1000 # mV
j=array([float(x) for x in '0 0.590 1.438 3.507'.split(' ')])*1000/10**4 # mA/cm^2

eta=e[1:]-e[0]
ln_j=log(j[1:])

h=array([ones(len(eta)),eta]).T # designmatrix

pseudoinverse=inv((h.T.dot(h))).dot(h.T)
params=pseudoinverse.dot(ln_j)

eta1=concatenate([[0],eta])
h1=concatenate([[[1,0]],h])

alpha=1-params[1]*1000*8.3145*293/96485

plt.figure()
plt.subplot(2,1,1)
plt.plot(eta,ln_j,'o',eta1,h1.dot(params))
plt.ylabel(r'$ln\left(j/(mA cm^{-2})\right)$')
plt.xlabel(r'$\eta / mV$')
plt.text(30,-3,
    '$j_0='+'{:0.5}'.format(exp(params[0]))+r'~mA\cdot cm^{-2}'+'$'+'\n'
    r'$\frac{(1-\alpha)F n}{R T}='+
    '{:0.5}'.format(params[1]*1000)+r'\frac{1}{V}\Rightarrow\alpha='+'{:0.5}'.format(alpha)+
    '$',ha='left')
plt.title('Atkins 8. E. Prob. 25.20')

plt.subplot(2,1,2)
eta_compl=linspace(-300,240,40)/1000 # V
j0=exp(params[0])
j_compl=j0*(exp((1-alpha)*96485*1*eta_compl/(8.3145*293))-exp(-alpha*96485*1*eta_compl/(8.3145*293)))
low_eta=array([-0.12,0.12])
j_low_eta=j0*low_eta*95485*1/(8.3145*293)
plt.plot(eta/1000,j[1:],'o',eta_compl,j_compl,'-',low_eta,j_low_eta,'--')
plt.xlabel(r'$\eta / V$')
plt.ylabel('$j/(mA cm^{-2})$')
plt.legend(['exp','Butler-Volmer',r'$-0,12~V\leq\eta\leq 0,12~V$'])
plt.tight_layout()

e=array([float(x) for x in '1.50 1.58 1.63 1.72 1.87 1.98 2.10'.split(' ')])*-1*1000 # mV
j=array([float(x) for x in '10 30 50 100 200 250 290'.split(' ')])*1000/10**4*-1 # mA/cm^2

e0=e[0] # overpotential at j~0 in mV
eta=e[2:]-e0 # 
ln_j=log(-j[2:])

h=array([ones(len(eta)),eta]).T # designmatrix

pseudoinverse=inv((h.T.dot(h))).dot(h.T)
params=pseudoinverse.dot(ln_j)

eta1=concatenate([[0],eta])
h1=concatenate([[[1,0]],h])

err=ln_j-h.dot(params)
r2=1-err.dot(err)/((ln_j-mean(ln_j)).dot(ln_j-mean(ln_j)))
rmsd=sqrt(err.dot(err)/len(err))

alpha=-params[1]*1000*8.3145*293/96485

plt.figure()
plt.subplot(2,2,1)
plt.plot(eta,ln_j,'o',eta1,h1.dot(params))
plt.ylabel(r'$ln\left(j/(mA cm^{-2})\right)$')
plt.xlabel(r'$\eta / mV$')
plt.text(-600,2,
    '$j_0='+'{:0.5}'.format(exp(params[0]))+r'~mA\cdot cm^{-2}'+'$'+'\n'
    r'$\frac{-\alpha \cdot F n}{R T}='+
    '{:0.5}'.format(params[1]*1000)+r'\frac{1}{V}\Rightarrow\alpha='+'{:0.5}'.format(alpha)+'$'+'\n'+
    '$'+'R^2='+'{:0.5}'.format(r2)+', RMSD='+'{:0.5}'.format(rmsd)+'$'
    ,ha='left')
plt.title('Atkins 8. E. Prob. 25.21')

eta_compl=linspace(-600,50,40)/1000 # V
j0=exp(params[0])
j_compl=j0*(exp((1-alpha)*96485*1*eta_compl/(8.3145*293))-exp(-alpha*96485*1*eta_compl/(8.3145*293)))
low_eta=array([-0.12,0.12])
j_low_eta=j0*low_eta*95485*1/(8.3145*293)
ax=plt.subplot(2,2,3)
plt.plot((e-e0)/1000,j,'o',eta_compl,j_compl,'-',low_eta,j_low_eta,'--')
plt.xlabel(r'$\eta / V$')
plt.ylabel('$j/(mA cm^{-2})$')
plt.legend(['exp','Butler-Volmer',r'$-0,12~V\leq\eta\leq 0,12~V$'],bbox_to_anchor=(0,1.02,1,.2),loc='lower left')
xticks=ax.get_xticks()
xlim=ax.get_xlim()
xlabels=xticks+e0/1000
secax=ax.twiny()
secax.set_xticks(xticks)
secax.set_xlim(xlim)
secax.set_xticklabels(xlabels)
secax.xaxis.set_ticks_position('bottom')
secax.xaxis.set_label_position('bottom')
secax.spines['bottom'].set_position(('outward',40))
secax.set_xlabel('E / V')


e0=e[0] # overpotential at j~0 in mV
eta=e[2:-1]-e0 # 
ln_j=log(-j[2:-1])

h=array([ones(len(eta)),eta]).T # designmatrix

pseudoinverse=inv((h.T.dot(h))).dot(h.T)
params=pseudoinverse.dot(ln_j)

eta1=concatenate([[0],eta])
h1=concatenate([[[1,0]],h])

err=ln_j-h.dot(params)
r2=1-err.dot(err)/((ln_j-mean(ln_j)).dot(ln_j-mean(ln_j)))
rmsd=sqrt(err.dot(err)/len(err))

alpha=-params[1]*1000*8.3145*293/96485

plt.subplot(2,2,2)
plt.plot(eta,ln_j,'o',eta1,h1.dot(params))
plt.ylabel(r'$ln\left(j/(mA cm^{-2})\right)$')
plt.xlabel(r'$\eta / mV$')
plt.text(-450,2,
    '$j_0='+'{:0.5}'.format(exp(params[0]))+r'~mA\cdot cm^{-2}'+'$'+'\n'
    r'$\frac{-\alpha \cdot F n}{R T}='+
    '{:0.5}'.format(params[1]*1000)+r'\frac{1}{V}\Rightarrow\alpha='+'{:0.5}'.format(alpha)+'$'+'\n'+
    '$'+'R^2='+'{:0.5}'.format(r2)+', RMSD='+'{:0.5}'.format(rmsd)+'$'
    ,ha='left')
plt.title('Atkins 8. E. Prob. 25.21')

eta_compl=linspace(-600,50,40)/1000 # V
j0=exp(params[0])
j_compl=j0*(exp((1-alpha)*96485*1*eta_compl/(8.3145*293))-exp(-alpha*96485*1*eta_compl/(8.3145*293)))
low_eta=array([-0.12,0.12])
j_low_eta=j0*low_eta*95485*1/(8.3145*293)
ax=plt.subplot(2,2,4)
plt.plot((e-e0)/1000,j,'o',eta_compl,j_compl,'-',low_eta,j_low_eta,'--')
plt.xlabel(r'$\eta / V$')
plt.ylabel('$j/(mA cm^{-2})$')
plt.legend(['exp','Butler-Volmer',r'$-0,12~V\leq\eta\leq 0,12~V$'],bbox_to_anchor=(0,1.02,1,.2),loc='lower left')
#secax=ax.secondary_xaxis('top',functions=(lambda x:x+e0/1000, lambda x:x-e0/1000))
#secax.set_xlabel('E / V')
xticks=ax.get_xticks()
xlim=ax.get_xlim()
xlabels=xticks+e0/1000
secax=ax.twiny()
secax.set_xticks(xticks)
secax.set_xlim(xlim)
secax.set_xticklabels(xlabels)
secax.xaxis.set_ticks_position('bottom')
secax.xaxis.set_label_position('bottom')
secax.spines['bottom'].set_position(('outward',40))
secax.set_xlabel('E / V')

plt.tight_layout()

plt.show()