## This is the user-script for plotting using gleaner tool.

# Path to gleaner (Better use environment variable PYTHONPATH!)
glrPath = '/home/gk779/apes/gleaner'
import sys
import os

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

## Import all required modules
import matplotlib.ticker as mtick
import matplotlib.pyplot as mplt
sys.path.append(glrPath)
import gleaner
import logging
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)

# Analytical solution
length = 0.5e-6
ref_pot = 10e-3
permit = 6.95e-10
moleDens0 = 1e-1 # 1 M = 1 mol/liter = 1e3 mol/m^3
valence_sqr = 1
charge = 1.6021e-19
faraday = 96485.3365
gasConst = 8.3144621
k_b = 1.3805e-23
temp = 273
N_A = 6.02e23 
from numpy import exp 
def analy_pot(x, moleDens):
  # reciprocal of debye length
  k = math.sqrt(2*moleDens*valence_sqr*faraday**2/(permit*gasConst*temp))*length
  #print('k ', k)
  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 

## -------------------------------------------------------------------------- ##
logging.basicConfig(level=logging.INFO)
## -------------------------------------------------------------------------- ##
logging.info('Started creating plots ...') 
 
import glob
import re
from operator import itemgetter

cases = ['C1','C01']
moleDensList = [1,0.1]
labels = ['$c_{\infty}=1$','$c_{infty}=0.1$']
import sqlite3
# data base filename
dbname = 'track.db'
# load database if exist else load tracking files and add to database
if os.path.isfile(dbname):
  print ('Processing data from existing database')
#  os.remove(dbname)
  sqlcon = sqlite3.connect(dbname)
else:
  print ('Processing data from tracking files')

  for iCase in cases:
    fname = glob.glob(iCase+'/*line_p00000_*.res')
    print ('Processing time file: *'+fname[0])

    # Load text, dump into a database with specific tabname to get columns later
    sqlcon = gleaner.tracking_to_db(fname = fname, \
                                      dbname=dbname, tabname='pot'+iCase)
## -------------------------------------------------------------------------- ##
print ('Potential over height:')
fig = mplt.figure()
ax = fig.add_subplot(111)
get_data_for_cols = ['coordX','potential_phy']

import numpy as np
x1 = np.linspace(0, length,50)
pot = analy_pot(x1, 0.1)
mplt.plot(x1/length, pot*1e3, 'x', color = 'r', label = 'Analytic, $c_{\infty}=0.1\ mol/m^3$')

# Simulation result
[x, y] = gleaner.get_columns(sqlcon, tabname='pot'+cases[0], \
                             columns=get_data_for_cols)
x, y = zip(*sorted(zip(x,y))) # sort of needed
xnew = []
ynew = []
for i in range(len(x)):
  xnew.append(x[i][0]/length)
  ynew.append(y[i][0]*1e3)
mplt.plot(xnew, ynew, '-', color='k',label = 'Simulation, $c_{\infty}=0.1\ mol/m^3$')

pot = analy_pot(x1, 1.0)
mplt.plot(x1/length, pot*1e3, 'o', color = 'r', label = 'Analytic, $c_{\infty}=1\ mol/m^3$')

# Simulation result
[x, y] = gleaner.get_columns(sqlcon, tabname='pot'+cases[1], \
                             columns=get_data_for_cols)
x, y = zip(*sorted(zip(x,y))) # sort of needed
xnew = []
ynew = []
for i in range(len(x)):
  xnew.append(x[i][0]/length)
  ynew.append(y[i][0]*1e3)
mplt.plot(xnew, ynew, '--', color='k',label = 'Simulation, $c_{\infty}=1\ mol/m^3$')

# plot setting
mplt.legend(loc=9, ncol=1,borderaxespad=0, \
            prop={'size':font_size}).get_frame().set_lw(0.0)
mplt.xlabel('y/H')
mplt.ylabel('Electric potential, $\psi$ ($kV$)')
mplt.grid(True,which="major",ls="-")
ax.yaxis.set_major_formatter(y_formatter)
mplt.xlim(0,1.0)

# save fig
figsize = [8,6]
fig = mplt.gcf()
fig.set_size_inches(figsize[0],figsize[1])
mplt.savefig('DoubleLayerPotentialOverHeight.png', dpi=100, format='png', \
             bbox_inches="tight",interpolation=None)
## -------------------------------------------------------------------------- ##
mplt.show()
logging.info('Plots created')
