plot  = []

## parameters
tau = 0.0001
mH2O = 18.01528e-3
mNa = 22.98977e-3
mCl = mNa
import math
xpos = 0.130900
diff_H2O_Na = 1.31570699e-9
diff_H2O_Cl = 2.097388e-9
omega = 0.95*math.pi
refNumDens = 1.0/1000.0
deltaRho1 = 0.001
deltaRho2 = 0.001
n0 = 1000.0
faraday = 96485.34
length = 2.0*math.pi
s = 2.0*math.pi/length
#Fickian diffusivity coeffcients
#D = 1.0 / ( (mNa/diff_H2O_Na)*(0.5-0.5*tau*deltaRho1-0.5*tau+deltaRho2) \
#  + (mNa/diff_H2O_Cl)*(0.5-0.5*tau*deltaRho1-0.5*tau-deltaRho2) )
D = 1.0 / (n0*((diff_H2O_Na+diff_H2O_Cl)/(2.0*diff_H2O_Na*diff_H2O_Cl)*(1.0 - deltaRho1*tau - tau) \
  + (diff_H2O_Cl-diff_H2O_Na)/(diff_H2O_Na*diff_H2O_Cl)*deltaRho2))
# constant of the velocity term
velConst = 2*deltaRho2*faraday*n0
rho0 = 22.98
# external electric field
E = 1e-3

show_plot = False
font_family= 'serif'
font_type= 'Times New Roman'
font_size = 16
figsize = [8,6]

infile = ['mix_point_*.res']
infile_bgk = ['tracking_bgk/mix_point_*.res']
infile_mrt = ['tracking_mrt/mix_point_*.res']

plot.append(dict(kind='xy', data=infile_bgk,
                 col=[1,2], npoints = 100, startplot = True,
                 label = 'LBM BGK', ls = 'k^', lw = 2))

plot.append(dict(kind='xy', data=infile_mrt,
                 col=[1,2], npoints = 90,
                 label = 'LBM MRT', ls = 'kd', lw = 2))

plot.append(dict(kind='function', data=infile_bgk,
                 col=[1,2,5], points = 50, endplot = True, 
                 fun = lambda *x: (math.exp((-1.0)*(s**2.0)*D*x[0])*
                 math.sin( s*xpos - s*(velConst/(2.0*x[2]))*E*x[0]**2.0 ) 
                 + (deltaRho1+1))*tau*n0,
                 xmin = 0, xmax = 1.5, ymin = 0.0, ymax = 0.25,
                 xlabel = 'time [s]', ylabel = 'concentration [$mol/m^3$]',
                 legendsize = 14, 
                 label = 'Analytic', ls = 'k', lw = 2))

plot.append(dict(kind='error', data=infile_bgk,
                 col=[1,2,5], points = 50, endplot = True, 
                 subplot = True, subaxes = [0.25,0.5,0.3,0.3],
                 reference = dict(kind='fun', 
                 fun = lambda *x: (math.exp((-1.0)*(s**2.0)*D*x[0])*
                 math.sin( s*xpos - s*(velConst/(2.0*x[2]))*E*x[0]**2.0 ) 
                 + (deltaRho1+1))*tau*n0),
                 label = 'BGK', ls = 'k-.'))

plot.append(dict(kind='error', data=infile_mrt,
                 col=[1,2,5], points = 50, endplot = True, 
                 subplot = True, subaxes = [0.25,0.5,0.3,0.3],
                 reference = dict(kind='fun', 
                 fun = lambda *x: (math.exp((-1.0)*(s**2.0)*D*x[0])*
                 math.sin( s*xpos - s*(velConst/(2.0*x[2]))*E*x[0]**2.0 ) 
                 + (deltaRho1+1))*tau*n0),
                 xmin = 0, xmax = 1.5, ymin = -0.005, ymax = 0.005,
                 #xlabel = 'time [$10^17s$]', ylabel = 'concentration [$mol/m^3$]',
                 #legendsize = 14, legend_loc = 1,
                 xticks = [0,0.5,1.0,1.5], yticks=[-0.005,-0.0025,0,0.0025,0.005],
                 title = 'Deviation', label = 'MRT',legend_loc=2, 
                 ls = 'k--', 
                 figname = 'dilute_analyticNumericComparisonFixedPos', format = 'png'
                 ))

