dia_inner = 2.77e-3 --m
dia_outer = 203.2e-3 --m
rad_inner = dia_inner/2.0
rad_outer = dia_outer/2.0
nDia_inner = 10
dx = dia_inner/nDia_inner
nDia_outer = math.ceil(dia_outer/dx)
nLength_bnd = 2*nDia_outer + 4
level = math.ceil(math.log(nLength_bnd)/math.log(2))
length_bnd = (2^level)*dx
--length_bnd = 2*dia_outer
--level = 9
--dx = length_bnd/2^level
--nLength = length_bnd/dx
qVal = true

dx_half = dx*0.5

dt =dx*dx
--print('dx ',dx)
--print('dt ',dt)
--print('nOD ', nDia_outer)
--print('nID ', nDia_inner)
--print('nOD-nID ', nDia_outer-nDia_inner)
permit = 8.854187e-12 
moleDens0 = 1e-3
valence_sqr = 1
charge = 1.60217657e-19
k_b = 1.3805e-23
temp = 273
N_A = 6.02e23 

ref_pot = 0.0 
pot_inner = 50e3 --V
pot_outer = 0.0
Q = 0--20e-6 -- C/m^3
-- Parameter for analytical solution
fac = Q/(4.0*permit)
function analy_pot(x,y,z,t)  
  r = math.sqrt(x^2+y^2+z^2)
  a = rad_inner
  b = rad_outer
  log_fac = (math.log(r)-math.log(b))/(math.log(a)-math.log(b))
  return fac*(b^2-r^2) + (pot_inner-fac*(b^2-a^2))*log_fac
end


