# Do not use X-server to create and save plot
import matplotlib
#matplotlib.use('Agg')

import numpy as np
import matplotlib.pyplot as mplt
import matplotlib.ticker as mtick
from numpy import exp 
import math

# font setting
from matplotlib import rc
font_size = 12
font_family = 'serif'
font_type = 'Times New Roman'
rc('text',usetex=True)
font = {'family':font_family,'%s'%font_family:font_type,'size':font_size}
rc('font',**font)

#axis without scientific notation
y_formatter = mtick.ScalarFormatter(useOffset=False)

length = 0.5e-6
nLength = 500
ref_pot = -10e-3
permit = 6.95e-10
moleDens0 = 1e-4*1e3 # 1 M = 1 mol/liter = 1e3 mol/m^3
valence_sqr = 1
charge = 1.6021e-19
k_b = 1.3805e-23
temp = 273
N_A = 6.02e23 
# reciprocal of debye length
lambda_inv = math.sqrt(2*moleDens0*N_A*valence_sqr*charge**2/(permit*k_b*temp))
k = lambda_inv*length
print(length/nLength, k, 1.0/lambda_inv)

def analytical_1(x):
  ek = exp(k)
  eminusk = exp(-k)
  term_1 = (ek - 1)/(ek-eminusk)*exp(-k*x/length)  
  term_2 = (1-eminusk)/(ek-eminusk)*exp(k*x/length)
  return (term_1+term_2)*ref_pot 

def analytical_2(x):
  fac_tanh = math.tanh(valence_sqr*charge*ref_pot/(4*k_b*temp))
  fac = 2.0*k_b*temp/(charge*valence_sqr)
  exp_lambda = exp(-x*lambda_inv)
  return fac*np.log( (1.0+fac_tanh*exp_lambda) / (1.0-fac_tanh*exp_lambda) )

def analytical_3(x):
  fac = valence_sqr*charge*ref_pot/(k_b*temp)
  exp_lambda = exp(-x*lambda_inv)
  #print(fac)
  #print(exp_lambda)
  return 2.0*np.log( (1.0+fac*exp_lambda) / (1.0-fac*exp_lambda) )

x = np.linspace(0,length,nLength)

fig = mplt.figure()
ax = fig.add_subplot(111)

analy_1 = analytical_1(x)
mplt.plot(x/length, analy_1, '-', color = 'r', label= 'Analy1')

analy_2 = analytical_2(x)
mplt.plot(x/length, analy_2, '-', color = 'b', label= 'Analy2')

#analy_3 = analytical_3(x)
#mplt.plot(x/length, analy_3, '-', color = 'g', label= 'Analy3')


# plot setting
mplt.legend(loc=3, ncol=1,borderaxespad=0, \
            prop={'size':font_size}).get_frame().set_lw(0.0)
mplt.xlabel('$x/Len$')
mplt.ylabel('Potential (V)')
mplt.grid(True,which="major",ls="-")
ax.yaxis.set_major_formatter(y_formatter)

# save fig
figsize = [8,6]
fig = mplt.gcf()
fig.set_size_inches(figsize[0],figsize[1])
mplt.savefig('PotentialOverX.png', dpi=100, format='png', \
             bbox_inches="tight",interpolation=None)

mplt.show()
